Solving a system of two coupled ordinary differential equations up to fourth-order

Question

I have the following (and terrible) system of ordinary differential equations up to fourth order (one of the equations is like a restriction for the other) for the one variable real functions $a(t)$ and $b(t)$:

$0=\frac{-12 \gamma (3 \alpha -\beta ) a(t)^2 b(t)^2 a''(t)^2-24 \gamma (3 \alpha -\beta ) a(t)^2 b(t) a'(t)^2 b''(t)+60 \gamma (3 \alpha -\beta ) a(t)^2 a'(t)^2 b'(t)^2-36 \gamma (3 \alpha -\beta ) b(t)^2 a'(t)^4+6 \gamma a(t)^2 b(t)^4 a'(t)^2+24 \gamma (3 \alpha -\beta ) a^{(3)}(t) a(t)^2 b(t)^2 a'(t)-24 b'(t) \left(\gamma (3 \alpha -\beta ) a(t) b(t) a'(t)^3+2 \gamma (3 \alpha -\beta ) a(t)^2 b(t) a'(t) a''(t)\right)+24 \gamma (3 \alpha -\beta ) a(t) b(t)^2 a'(t)^2 a''(t)-\Lambda a(t)^4 b(t)^2}{2 \gamma a(t)^4 b(t)^4}$

$0=-\frac{8 \gamma (3 \alpha -\beta ) a(t)^3 a^{(4)}(t) b(t)^3+12 \gamma (3 \alpha -\beta ) a(t)^2 b(t)^3 a''(t)^2-8 \gamma (3 \alpha -\beta ) a(t)^3 b^{(3)}(t) b(t)^2 a'(t)-120 \gamma (3 \alpha -\beta ) a(t)^3 a'(t) b'(t)^3+12 \gamma (3 \alpha -\beta ) b(t)^3 a'(t)^4+2 \gamma a(t)^2 b(t)^5 a'(t)^2+16 \gamma (3 \alpha -\beta ) a(t)^2 a^{(3)}(t) b(t)^3 a'(t)+60 b'(t)^2 \left(2 \gamma (3 \alpha -\beta ) a(t)^3 b(t) a''(t)+\gamma (3 \alpha -\beta ) a(t)^2 b(t) a'(t)^2\right)-16 b''(t) \left(2 \gamma (3 \alpha -\beta ) a(t)^3 b(t)^2 a''(t)-5 \gamma (3 \alpha -\beta ) a(t)^3 b(t) a'(t) b'(t)+\gamma (3 \alpha -\beta ) a(t)^2 b(t)^2 a'(t)^2\right)+4 a''(t) \left(\gamma a(t)^3 b(t)^5-12 \gamma (3 \alpha -\beta ) a(t) b(t)^3 a'(t)^2\right)-4 b'(t) \left(12 \gamma (3 \alpha -\beta ) a(t)^3 a^{(3)}(t) b(t)^2-12 \gamma (3 \alpha -\beta ) a(t) b(t)^2 a'(t)^3+\gamma a(t)^3 b(t)^4 a'(t)+18 \gamma (3 \alpha -\beta ) a(t)^2 b(t)^2 a'(t) a''(t)\right)-\Lambda b(t)^7}{2 \gamma a(t)^2 b(t)^7}$

Where I am using prime to denote derivative with respect to time $t$. In Mathematica they look like this

Ec1 = 1/(
  2 \[Gamma] a[t]^4 b[t]^4) (6 \[Gamma]*a[t]^2 b[t]^4 a'[t]^2 - 
    36 (3 \[Alpha] - \[Beta]) \[Gamma]*b[t]^2 a'[t]^4 + 
    24 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t] b[t]^2 a'[t]^2 a''[t] - 
    12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]^2 a''[t]^2 + 
    24 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]^2 a'[t] a'''[t] + 
    60 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 a'[t]^2 b'[t]^2 - 
    24 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t] a'[t]^2 b''[t] - 
    a[t]^4 b[t]^2 \[CapitalLambda] - 
    24 ((3 \[Alpha] - \[Beta]) \[Gamma]*a[t] b[t] a'[t]^3 + 
       2 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t] a'[t] a''[t]) b'[
      t])

Ec2 = -(1/(
   2 \[Gamma] a[t]^2 b[t]^7)) (2 \[Gamma]*a[t]^2 b[t]^5 a'[t]^2 + 
    12 (3 \[Alpha] - \[Beta]) \[Gamma]*b[t]^3 a'[t]^4 + 
    12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]^3 a''[t]^2 + 
    16 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2*b[t]^3 a'[t]*a'''[t] + 
    8 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^3 a''''[t] - 
    120 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3 a'[t]*b'[t]^3 - 
    8 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^2*a'[t]*b'''[t] + 
    b[t]^7 (-\[CapitalLambda]) + 
    60 ((3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]*a'[t]^2 + 
       2 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3 b[t]*a''[t]) b'[
       t]^2 + 4 (\[Gamma]*a[t]^3 b[t]^5 - 
       12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]*b[t]^3*a'[t]^2) a''[
      t] - 4 (\[Gamma]*a[t]^3 b[t]^4 a'[t] - 
       12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]*b[t]^2*a'[t]^3 + 
       18 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2*
        b[t]^2 a'[t] a''[t] + 
       12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^2 a'''[t]) b'[
      t] - 16 ((3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2*b[t]^2*
        a'[t]^2 + 
       2 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^2*a''[t] - 
       5 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]*a'[t]*b'[t]) b''[
      t])

[Alpha] = 1/20; [Beta] = 1/5; [Gamma] = 1; [CapitalLambda] = 1; 

sol = NDSolve[{Ec1 == 0,Ec2== 0, ic}, {a[t], b[t]}, {t, 0,1} ] 

Here, $\alpha$,$\beta$,$\gamma$ and $\Lambda$ are constants of $O(1)$. Of course, a lot of initial conditions need to be given. Some reasonable initial conditions from the context of the problem are $a(t=1)=b(t=1)=1$, and their first derivatives should be positive $a'(t=1),b'(t=1)>0$ but not so big (0.5, for example). I don't have more information in order to give solid initial conditions for the higher derivatives, but it is reasonable again to choose them positive and the initial condition for the second derivative $b''(t=1)$ can be obtained from the first equation by consistency. One could choose, for example:

ic = {a[1] == 1, a'[1] == 1, a''[1] == 1/2, a'''[1]==1/2 , b[1] == 1, 
 b'[1] == 1, b''[1] == -(103/24)}

The problem is that Mathematica is cannot solve this system, giving errors like "NDSolve::ndcf: Repeated convergence test failure at t == 1.`; unable to continue.".

I have tried to use AsymptoticDSolveValue too (https://reference.wolfram.com/language/ref/AsymptoticDSolveValue.html), but it doesn't work. I would like to know if it's possible to know something about the solution of this system, if one can check its singularities, behaviour near $t=0$, etc.

Answer 1

This problem can be solved as optimization problem with using the Euler wavelets collocation method as follows

eq1 = 1/(2 \[Gamma] a[t]^4 b[t]^4) (6 \[Gamma]*
      a[t]^2 b[t]^4 a'[t]^2 - 
     36 (3 \[Alpha] - \[Beta]) \[Gamma]*b[t]^2 a'[t]^4 + 
     24 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t] b[t]^2 a'[t]^2 a''[t] - 
     12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]^2 a''[t]^2 + 
     24 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]^2 a'[t] a'''[t] + 
     60 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 a'[t]^2 b'[t]^2 - 
     24 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t] a'[t]^2 b''[t] - 
     a[t]^4 b[t]^2 \[CapitalLambda] - 
     24 ((3 \[Alpha] - \[Beta]) \[Gamma]*a[t] b[t] a'[t]^3 + 
        2 (3 \[Alpha] - \[Beta]) \[Gamma]*
         a[t]^2 b[t] a'[t] a''[t]) b'[t]);

eq2 = -(1/(2 \[Gamma] a[t]^2 b[t]^7)) (2 \[Gamma]*
      a[t]^2 b[t]^5 a'[t]^2 + 
     12 (3 \[Alpha] - \[Beta]) \[Gamma]*b[t]^3 a'[t]^4 + 
     12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]^3 a''[t]^2 + 
     16 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2*b[t]^3 a'[t]*a'''[t] + 
     8 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^3 a''''[t] - 
     120 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3 a'[t]*b'[t]^3 - 
     8 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^2*a'[t]*b'''[t] + 
     b[t]^7 (-\[CapitalLambda]) + 
     60 ((3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2 b[t]*a'[t]^2 + 
        2 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3 b[t]*a''[t]) b'[
        t]^2 + 4 (\[Gamma]*a[t]^3 b[t]^5 - 
        12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]*b[t]^3*a'[t]^2) a''[
       t] - 4 (\[Gamma]*a[t]^3 b[t]^4 a'[t] - 
        12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]*b[t]^2*a'[t]^3 + 
        18 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2*
         b[t]^2 a'[t] a''[t] + 
        12 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^2 a'''[t]) b'[
       t] - 16 ((3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^2*b[t]^2*
         a'[t]^2 + 
        2 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]^2*a''[t] - 
        5 (3 \[Alpha] - \[Beta]) \[Gamma]*a[t]^3*b[t]*a'[t]*b'[t]) b''[
       t]);

\[Alpha] = 1/20; \[Beta] = 1/5; \[Gamma] = 1; \[CapitalLambda] = 1;

Let check that initial conditions ic = {a[1] == 1, a'[1] == 1, a''[1] == 1/2, a'''[1] == 1/2, b[1] == 1,b'[1] == 1, b''[1] == -515/120} satisfier to eq1, we have

 eq1 /. {a[t] -> 1, a'[t] -> 1, a''[t] -> 1/2, a'''[t] -> 1/2, 
  b[t] -> 1, b'[t] -> 1, b''[t] -> -515/120}

Out[]= 0 

Then we transform eq1, eq2 to the system of algebraic equations in a form

ic = {a[1] == 1, a'[1] == 1, a''[1] == 1/2, a'''[1] == 1/2, b[1] == 1,
    b'[1] == 1, b''[1] == -515/120};
eqs = {eq1 == 0, eq2 == 0};
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 = 4;
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]^4 b[t]^4 eq1/18 /. rule; eqn2 = 
 a[t]^2 b[t]^7 eq2/30 /. 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}]; 

Numerical solution

sol = NMinimize[{eqs1 . eqs1 + eqs2 . eqs2, bc}, var, 
  MaxIterations -> 1000, Method -> "DifferentialEvolution"]

Visualization solution a[t], b[t] with collocation points

Show[Plot[Evaluate[{a0[t], b0[t]} /. sol[[2]]], {t, 0, 1}, 
  PlotLegends -> {"a", "b"}, AxesLabel -> Automatic], 
 ListPlot[{Table[{t, a0[t]}, {t, tcol}], 
    Table[{t, b0[t]}, {t, tcol}]} /. sol[[2]], 
  PlotStyle -> {Blue, Red}]]

To compute solution at t>1 we use

ic1 = {a[0] == 1, a'[0] == 1, a''[0] == 1/2, a'''[0] == 1/2, 
  b[0] == 1, b'[0] == 1, b''[0] == -515/120}; rul0 = 
 rule /. t -> 0; bc1 = ic1 /. rul0;
sol1 = NMinimize[{eqs1 . eqs1 + eqs2 . eqs2, bc1}, var, 
  MaxIterations -> 1000, Method -> "DifferentialEvolution"]

Visualization

Show[Plot[Evaluate[{a0[t - 1], b0[t - 1]} /. sol1[[2]]], {t, 1, 2}, 
  PlotLegends -> {"a", "b"}, AxesLabel -> Automatic], 
 ListPlot[{Table[{t + 1, a0[t]}, {t, tcol}], 
    Table[{t + 1, b0[t]}, {t, tcol}]} /. sol1[[2]], 
  PlotStyle -> {Blue, Red}]]