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

NDSolve
needs a complete set of initial conditions, in your case 4*2. By the way it is good praxis toRationalize
the equations before solving $\endgroup$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 ineq2
, so I think you have a problem there in that the highest order derivatives are not determined by the lower order ones. MaybeNDSolve
can solve the system as a DAE, but I doubt it. Maybe there's a mistake in the set up? $\endgroup$Reduce[Rationalize@{eq1, D[eq2, t]}, Derivative[4][a][t]]
showsa
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}
. $\endgroup$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])
$\endgroup$