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

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.

• solve this system, giving errors like can you show the full code you used? Including the NDSolve command? Commented Jun 14, 2023 at 7:39
• Your code is wrong. You write \[CapitalLambda][t] in one place then write \[CapitalLambda] on its own with no argument in another place. Besides this, how is NDSolve supposed to solve this if it does not know what this function is? You also have \[Alpha] , \[Beta] in there with no numerical values. This will not work for numerical solver. Commented Jun 14, 2023 at 7:52
• The main problem is that NDSolve says Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations. May be someone else will have an idea why it says so for your system. Commented Jun 14, 2023 at 8:57
• You make it inconvenient and tedious to work with your code: Correct your code above regarding $\Lambda$, write all the code in one copy-and-paste block including setting up the initial conditions, include the full NDSolve code, and in a comment block, the error message you're getting so all someone has to do is copy and paste and run it. Normally good practice to Latex the DE but yours is so messy it's only a distraction. Consider removing it.
– josh
Commented Jun 14, 2023 at 12:23
• If you have a system of ODEs where one of the equations is a constraint on the lower derivatives, then it will go better for you if you can figure out what the true degrees of freedom are for the system before you feed it into a numerical solver. You may also need to delve into the detailed documentation on numerical ODE solvers in Mathematica, particularly the section on DAE solvers. Godspeed. Commented Jun 14, 2023 at 15:56

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


• That looks fantastic. Thank you so much. I have a last question, is there any way to check this is the solution to the system? I am thinking in something like eqs /. sol /. {t -> RandomReal[]} to correctly check that the system is solved succesfully. Commented Jun 17, 2023 at 7:17
• @Axionlikeparticles This is not solution in ordinary sense, this is optimal solution with an average error of Sqrt[sol[[1]]/(2 nn)] for t<=1, and Sqrt[sol1[[1]]/(2 nn)] for t>1. You can also check bc/.sol[[2]] and bc1/.sol1[[2]]`. Commented Jun 17, 2023 at 8:15
• Yeah, but that checks are purely for the initial conditions. I guess that you mean that this is not the solution in the ordinary sense but it is an accurate numerical solution for the system, right? Commented Jun 17, 2023 at 16:35
• @Axionlikeparticles This is solution for the system with minimal error as possible, but it is not unique. Commented Jun 18, 2023 at 1:49
• Well yeah, but in my mind if it's the most accurate numerical one then it is the most important. Even more, expanding in series for example (as the other problem we discussed) gives a similar answer, so it is kinda reliable as the true solution. Thank you! I'm wondering one last thing, is it possible to do a series expansion without going trough the Euler wavepackets method and obtaining a similar solution? Commented Jun 18, 2023 at 8:08