# Solving a system of two fourth-order ordinary differential equations

I have the following fourth order ordinary differential equations for the functions $$a(t)$$ and $$b(t)$$:

eq1 = b[t]*
Derivative[1][a][t]^2*(-3.*Derivative[1][a][t]*Derivative[1][b][t] +
3.*b[t]*Derivative[2][a][t]) +
a[t]*(1.*b[t]^5 + 0.5*b[t]^4*Derivative[1][a][t]^2 -
3.75*Derivative[1][a][t]^2*Derivative[1][b][t]^2 +
b[t]*Derivative[1][a][t]*(4.5*Derivative[1][b][t]*Derivative[2][a][t] +
1.*Derivative[1][a][t]*Derivative[2][b][t]) +
b[t]^2*(-0.75*Derivative[2][a][t]^2 -
1.*Derivative[1][a][t]*Derivative[3][a][t])) +
a[t]^2*(-1.*b[t]^3*Derivative[1][a][t]*
Derivative[1][b][t] + (7.5*Derivative[1][a][t]*Derivative[1][b][t]^3)/b[t] +
1.*b[t]^4*Derivative[2][a][t] +
Derivative[1][b][t]*(-7.5*Derivative[1][b][t]*Derivative[2][a][t] -
5.*Derivative[1][a][t]*Derivative[2][b][t]) +
b[t]*(2.*Derivative[2][a][t]*Derivative[2][b][t] +
3.*Derivative[1][b][t]*Derivative[3][a][t] +
0.5*Derivative[1][a][t]*Derivative[3][b][t]) -
0.5*b[t]^2*Derivative[4][a][t]) == (0.75*b[t]^2*
Derivative[1][a][t]^4)/a[t]

eq2 = ((18.*b[t]^2*Derivative[1][a][t]^4)/a[t] -
12.*b[t]*
Derivative[1][a][t]^2*(-1.*Derivative[1][a][t]*Derivative[1][b][t] +
b[t]*Derivative[2][a][t]) +
6*a[t]*(2*b[t]^4*Derivative[1][a][t]^2 -
5.*Derivative[1][a][t]^2*Derivative[1][b][t]^2 +
2.*b[t]*Derivative[1][a][
t]*(2*Derivative[1][b][t]*Derivative[2][a][t] +
Derivative[1][a][t]*Derivative[2][b][t]) +
1.*b[t]^2*(Derivative[2][a][t]^2 -
2.*Derivative[1][a][t]*Derivative[3][a][t])))/b[t]^6 == -24


In order to solve them, I use the following code

NDSolve[{eq1, eq2, a[1] == 1, b[1] == 1}, {a, b}, {t, 0.5, 10}]


Of course, I have tried adding more initial conditions and changed the time interval but none of this works. Mathematica can't solve it and tells me the following:

"NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations."

"NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions."

How can I tell Mathematica to solve this system?

• NDSolve needs a complete set of initial conditions, in your case 4*2. By the way it is good praxis to Rationalize the equations before solving Apr 27 at 10:38
• Try Solve[{eq1, eq2}, {Derivative[4][a][t], Derivative[3][b][t]}], which tries to solve for the highest order derivative that I can seen. Neither occurs in eq2, so I think you have a problem there in that the highest order derivatives are not determined by the lower order ones. Maybe NDSolve can solve the system as a DAE, but I doubt it. Maybe there's a mistake in the set up? Apr 27 at 13:25
• I finally could experiment: Reduce[Rationalize@{eq1, D[eq2, t]}, Derivative[4][a][t]] shows a needs to be a constant function. But that implies (eq2 /. {a -> Function[t, C[1]]}) the second equation is inconsistent. I think there's a coding mistake somewhere in {eq1, eq2}. Apr 28 at 12:20
• @AlexTrounev Kind of, it is an extension from general relativity, where instead of the Einstein-Hilbert lagrangian density $\mathcal L = \sqrt{-g}\gamma R$ one has the lagrangian density of higher curvature terms (up to fourth order) given by $\mathcal L = \sqrt{-g}(\gamma R+\alpha R_{\mu\nu}R^{\mu\nu}-\beta R^2$. If one inserts the FLRW metric $ds^2=-b^2(t)dt^2+a^2(t)(dx^2+dy^2+dz^2)$ and applies the Euler-Lagrange equations for higher derivative Lagrangians, one obtains both equations above for $a$ and $b$. That's the context, but I don't think it helps much for the problem stated. Apr 29 at 8:31
• @AlexTrounev Sure, sorry for the delay. The lagrangian as a function of $a$ and $b$ is Lagrangian = 6/(a[t]*b[t]^5)*(\[Gamma]*a[t]^2*b[t]^4*a'[t]^2 + 2*(\[Alpha] - 3*\[Beta])*b[t]^2*a'[t]^4 + 2*(\[Alpha] - 3*\[Beta])*a[t]^2*b[t]^2*a''[t]^2 + 2*(\[Alpha] - 3*\[Beta])*a[t]^2*a'[t]^2* b'[t]^2 + (\[Gamma]*a[t]^3*b[t]^4 + 2*(\[Alpha] - 6*\[Beta])*a[t]*b[t]^2*a'[t]^2)*a''[t] - (\[Gamma]*a[t]^3*b[t]^3*a'[t] + 2*(\[Alpha] - 6*\[Beta])*a[t]*b[t]*a'[t]^3 + 4*(\[Alpha] - 3*\[Beta])*a[t]^2*b[t]*a'[t]*a''[t])*b'[t]) Apr 30 at 9:15

As it mentioned by Axionlikeparticles, the system of equations has been derived from Lagrangian $$\mathcal L = \sqrt{-g}(\gamma R+\alpha R_{\mu\nu}R^{\mu\nu}-\beta R^2)$$ in some extension of GR in a case of the FLRW metric $$ds^2=-b^2(t)dt^2+a^2(t)(dx^2+dy^2+dz^2)$$. To reproduce this Lagrangian we use code

n = 4;
coord = {t, x, y, z};
metric = {{-b[t]^2, 0, 0, 0}, {0, a[t]^2, 0, 0}, {0, 0, a[t]^2,
0}, {0, 0, 0, a[t]^2}};

inversemetric = Simplify[Inverse[metric]];

affine :=
affine =
Simplify[
Table[(1/2)*
Sum[inversemetric[[i,
s]]*(D[metric[[s, j]], coord[[k]]] +
D[metric[[s, k]], coord[[j]]] -
D[metric[[j, k]], coord[[s]]]), {s, 1, n}], {i, 1, n}, {j,
1, n}, {k, 1, n}]];
riemann :=
riemann =
Simplify[
Table[D[affine[[i, j, l]], coord[[k]]] -
D[affine[[i, j, k]], coord[[l]]] +
Sum[affine[[s, j, l]]*affine[[i, k, s]] -
affine[[s, j, k]]*affine[[i, l, s]], {s, 1, n}], {i, 1,
n}, {j, 1, n}, {k, 1, n}, {l, 1, n}]];
ricci :=
ricci = Simplify[
Table[Sum[riemann[[i, j, i, l]], {i, 1, n}], {j, 1, n}, {l, 1, n}]]

scalar =
Simplify[
Sum[inversemetric[[i, j]]*ricci[[i, j]], {i, 1, n}, {j, 1, n}]];
Ril2 = Sum[
ricci[[i, k]] inversemetric[[i, l]] inversemetric[[k, m]] ricci[[l,
m]], {i, 1, 4}, {k, 1, 4}, {l, 1, 4}, {m, 1, 4}] // Simplify;

R2 = scalar^2;

L = Sqrt[-Det[
metric]] (\[Gamma] scalar + \[Alpha] Ril2 - \[Beta] R2) //
Simplify


As output we have

L=1/(a[t]^4 b[t]^6)
Sqrt[a[t]^6 b[
t]^2] (-\[Beta] (6 a[t] Derivative[1][a][t] Derivative[1][b][t] -
6 b[t] (Derivative[1][a][t]^2 +
a[t] (a^\[Prime]\[Prime])[t]))^2 + \[Gamma] a[t]^2 b[
t]^3 (-6 a[t] Derivative[1][a][t] Derivative[1][b][t] +
6 b[t] (Derivative[1][a][t]^2 +
a[t] (a^\[Prime]\[Prime])[t])) + \[Alpha] (9 a[
t]^2 (Derivative[1][a][t] Derivative[1][b][t] -
b[t] (a^\[Prime]\[Prime])[t])^2 +
3 (a[t] Derivative[1][a][t] Derivative[1][b][t] -
b[t] (2 Derivative[1][a][t]^2 +
a[t] (a^\[Prime]\[Prime])[t]))^2));


From this Lagrangian we can derive equations using $$\frac{\partial L}{\partial q} - \frac{d}{dt}\frac{\partial L}{\partial \dot q} + \frac{d^2}{dt^2}\frac{\partial L}{\partial \ddot q} = 0$$ for $$q=a(t)$$. Therefore we have

eq1 =
D[L, a[t]] - D[D[L, a'[t]], t] + D[D[L, a''[t]], t, t] // Simplify;

eq2 = D[L, b[t]] - D[D[L, b'[t]], t] // Simplify;


This system is not differ from that discussed above, so the problem of NDSolve usage still remains. For example,

eqs = {eq1 == 0, eq2 == 0}; ic = {a[1] == 1, a'[1] == 1,
a''[1] == 0.5, a'''[1] == 0.5, b[1] == 1, b'[1] == 1, b''[1] == 0.5};

sols = NDSolve[{eqs /. {\[Alpha] -> 1/2, \[Beta] -> 1/5, \[Gamma] ->
1}, ic}, {a, b}, {t, 0, 1}]


As output we have messages

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

NDSolve::parpiv: Zero pivot was detected during the numerical factorization or there was a problem in the iterative refinement process. It is possible that the matrix is ill-conditioned or singular.


The problem here is that eq2 is a constraint on eq1, not an independent equation. It means, that we can try to solve the system eq1,eq2,ic as an optimization problem. For this we need some ic consistent with eq2. Let check that

eq2 /. {a[t] -> 1, a'[t] -> 1, a''[t] -> 1/2, a'''[t] -> 1/2,
b[t] -> 1,
b'[t] -> 1, b''[t] -> -(21/8)} /. {\[Alpha] ->
1/2, \[Beta] -> 1/5, \[Gamma] -> 1}

(*Out[]= 0*)


Therefore our new initial conditions are given by

ic = {a[1] == 1, a'[1] == 1, a''[1] == 1/2, a'''[1] == 1/2, b[1] == 1,
b'[1] == 1, b''[1] == -21/8};


To convert eq1,eq2,ini into a system of algebraic equations we use the Euler wavelets collocation method described in our paper

OEm[m_, x_] :=
Sqrt[2 m +
1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0,
m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7;
With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dt = 1/(nn); tl = Table[l*dt, {l, 0, nn}]; tcol =
Table[(tl[[l - 1]] + tl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Int3 = Integrate[Int2, t1];
Int4 = Integrate[Int3, t1];
Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
int3[y_] := Int3 /. t1 -> y;
int4[y_] := Int4 /. t1 -> y;
A = Array[aa, nn]; B = Array[bb, nn];
a4[t_] := A . Psi[t];
a3[t_] := A . int1[t] + c1 ;
a2[t_] := A . int2[t] + c1 t + c2;
a1[t_] := A . int3[t] + c1 t^2/2 + c2 t + c3;
a0[t_] := A . int4[t] + c1 t^3/6 + c2 t^2/2 + c3 t + c4;
b3[t_] := B . Psi[t]; b2[t_] := B . int1[t] + d1 ;
b1[t_] := B . int2[t] + d1 t + d2;
b0[t_] := B . int3[t] + d1 t^2/2 + d2 t + d3;
rule = {Derivative[4][a][t] -> a4[t], Derivative[3][a][t] -> a3[t],
Derivative[2][a][t] -> a2[t], Derivative[1][a][t] -> a1[t],
a[t] -> a0[t], Derivative[4][b][t] -> b4[t],
Derivative[3][b][t] -> b3[t], Derivative[2][b][t] -> b2[t],
Derivative[1][b][t] -> b1[t], b[t] -> b0[t]};
eqn1 = a[t]^2 b[t]^6 eq1/6 /. rule; eqn2 = a[t] b[t]^6 eq2/6 /. rule;

rul1 = rule /. t -> 1; bc = ic /. rul1; eqs1 =
Flatten[Table[eqn1 /. t -> tcol[[i]], {i, nn}]]; eqs2 =
Flatten[Table[eqn2 /. t -> tcol[[i]], {i, nn}]]; var =
Join[A, B, {c1, c2, c3, c4, d1, d2, d3}];


The solution of optimization problem is given by

sol = NMinimize[{eqs1 . eqs1 + eqs2 . eqs2,
bc} /. {\[Alpha] -> 1/2, \[Beta] -> 1/5, \[Gamma] -> 1}, var,
MaxIterations -> 1000, Method -> "DifferentialEvolution"]

(*Out[]= {6.0224*10^-7, {aa[1] -> 0.662598, aa[2] -> 6.6182,
aa[3] -> 0.411078, aa[4] -> -4.25886, aa[5] -> 6.38316,
aa[6] -> 3.1353, aa[7] -> 3.66257, aa[8] -> -3.99448,
aa[9] -> 3.72062, aa[10] -> 1.33545, aa[11] -> 17.4548,
aa[12] -> 3.23394, aa[13] -> -2.96937, aa[14] -> 0.864018,
aa[15] -> -1.77324, aa[16] -> 1.55967, aa[17] -> -4.07709,
aa[18] -> -0.890433, aa[19] -> -1.54841, aa[20] -> -4.37257,
aa[21] -> 0.199925, aa[22] -> 6.07787, aa[23] -> 0.84347,
aa[24] -> -6.724, aa[25] -> 3.21171, aa[26] -> -0.85385,
aa[27] -> 0.273778, aa[28] -> 0.0249634, bb[1] -> -1.18045,
bb[2] -> 5.51932, bb[3] -> 0.613343, bb[4] -> -0.346739,
bb[5] -> 2.66884, bb[6] -> -10.0235, bb[7] -> -1.77049,
bb[8] -> 1.2248, bb[9] -> 0.599704, bb[10] -> 0.700833,
bb[11] -> 9.03793, bb[12] -> 2.85301, bb[13] -> -1.2703,
bb[14] -> 0.434722, bb[15] -> -4.6445, bb[16] -> -1.74593,
bb[17] -> -2.82947, bb[18] -> -1.06695, bb[19] -> -1.26808,
bb[20] -> -3.13756, bb[21] -> -0.275597, bb[22] -> -4.31528,
bb[23] -> 1.52497, bb[24] -> -5.6961, bb[25] -> 2.06735,
bb[26] -> -0.336646, bb[27] -> 0.161848, bb[28] -> 0.0455251,
c1 -> 0.013626, c2 -> 1.85716, c3 -> -0.235057, c4 -> 0.481964,
d1 -> 1.83272, d2 -> 0.654357, d3 -> -0.0931884}}*)


The mean error of optimal solution is about 10^-4. It is not perfect but we can check that

bc /. sol[[2]]

Out[]= {True, True, True, True, True, True, True}


Visualization

plot1=Plot[Evaluate[{a0[t], b0[t]} /. sol[[2]]], {t, 0, 1},
PlotLegends -> {"a", "b"}, AxesLabel -> Automatic]


Using modified ic we can extend this solution for t>1 as follows

ic1 = {a[0] == 1, a'[0] == 1, a''[0] == 1/2, a'''[0] == 1/2,
b[0] == 1, b'[0] == 1, b''[0] == -21/8};
bc = ic1 /. rul1;


This solution looks like

Update 1. We also can use series expansions in a form $$a=\sum_{k=0}^n a_kt^k, b=\sum_{k=0}^n b_kt^k$$, for example

 L = 6/(a[t]*b[t]^5)*(\[Gamma]*a[t]^2*b[t]^4*a'[t]^2 +
2*(\[Alpha] - 3*\[Beta])*b[t]^2*a'[t]^4 +
2*(\[Alpha] - 3*\[Beta])*a[t]^2*b[t]^2*a''[t]^2 +
2*(\[Alpha] - 3*\[Beta])*a[t]^2*a'[t]^2*
b'[t]^2 + (\[Gamma]*a[t]^3*b[t]^4 +
2*(\[Alpha] - 6*\[Beta])*a[t]*b[t]^2*a'[t]^2)*
a''[t] - (\[Gamma]*a[t]^3*b[t]^3*a'[t] +
2*(\[Alpha] - 6*\[Beta])*a[t]*b[t]*a'[t]^3 +
4*(\[Alpha] - 3*\[Beta])*a[t]^2*b[t]*a'[t]*a''[t])*b'[t]);

eq1 = D[L, a[t]] - D[D[L, a'[t]], t] + D[D[L, a''[t]], t, t] //
Simplify

eq2 = D[L, b[t]] - D[D[L, b'[t]], t] + D[D[L, b''[t]], t, t] //
Simplify

ic = {a[1] == 1, a'[1] == 1, a''[1] == 1/2, a'''[1] == 1/2, b[1] == 1,
b'[1] == 1, b''[1] == -21/8};
n = 11; A = Array[aa, n]; B = Array[bb, n]; T =
Table[(t - 1)^k, {k, 0, n - 1}]; a[t_] = A . T;
b[t_] = B . T; grid = Range[2/5, 1.4, 1/10];
eqs1 = Table[a[t] b[t]^6/6 eq1, {t, grid}]; eqs2 =
Table[a[t] b[t]^6/6 eq2, {t, grid}];


The optimal solution is given by

sol1 =
NMinimize[{eqs1 . eqs1 + eqs2 . eqs2,
ic} /. {\[Alpha] -> 1/2, \[Beta] -> 1/5, \[Gamma] -> 1},
Join[A, B], MaxIterations -> 1000, Method -> "DifferentialEvolution"]

(*Out[]= {0.0000233149, {aa[1] -> 1., aa[2] -> 1., aa[3] -> 0.25,
aa[4] -> 0.0833333, aa[5] -> 0.438961, aa[6] -> -0.30387,
aa[7] -> -0.317073, aa[8] -> 0.125383, aa[9] -> -0.182751,
aa[10] -> 0.374697, aa[11] -> 0.0144105, bb[1] -> 1., bb[2] -> 1.,
bb[3] -> -1.3125, bb[4] -> -1.13918, bb[5] -> 0.613892,
bb[6] -> 0.413397, bb[7] -> 0.466893, bb[8] -> -0.199954,
bb[9] -> 0.622894, bb[10] -> -0.567313, bb[11] -> 0.0127384}}*)


We can compare this solution (doted lines) to solution above computed with Euler wavelets (solid lines) for $$t\le 1$$ (for $$t>1$$ the correlation is not so good)

Show[plot1,
Plot[Evaluate[{a[t], b[t]} /. sol1[[2]]], {t, 0.4, 1.},
PlotStyle -> {{Blue, Dashed}, {Red, Dashed}}]]


• Thank you so much for your answer. I don't understand the algorithm very well, but it seems (you made a nice debugging there) you obtained a valid solution. Is there any way to avoid using the Euler wavelets collocation method? It seems it works fine, but it is tedious if one wants to set up new initial conditions several times. May 2 at 15:56
• @Axionlikeparticles If solution looks nice, then I can organize code in one function with parameters and region specification. May 2 at 17:39
• Yeah, the solution above is physically pathological because one change of signature, but that is because there's an ==0 in the equations and not an ==something (which has to do with energy content), which is very easy to change. I have tried it in another realistic case and it looks aparently nice. Your proposal is fantastic, I really appreciate it, and that would be awesome. Also, I have a question, why is the method used only valid between 0 and 1 and the ic need to be changed in order to obtain t>1? May 2 at 19:28
• @Axionlikeparticles The reason is the numerical solution convergence. The Euler wavelets are normalized and orthogonalized on the unit interval $0<=t<=1$. Also for $0<=t<=1$ we have better result compare to $1<t<=2$. It is why we not mixed two solutions. In practice we can map any interval $0<=t<=T$ to unit interval. Please note that your problem is autonomous, hence we can map $1<t<=2$ to the unit interval with substitution $t\rightarrow t-1$, it is why new ic stated at t=0. May 3 at 2:16
• @Axionlikeparticles See Update 1 to my answer. May 4 at 4:05