Solving a system of two fourth-order ordinary differential equations

Question

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?

Answer 1

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}}]]