Highly coupled nonlinear second order differential equations

Question

I have an action given by, $$S = \int dx \frac{1}{z^d} \sqrt{-f(z,u) u'^2 - 2 u' z' +1}$$

The dependent variables are $u$ and $z$ for which they are dependent on the parameter $x$. The equation of motion (EOM) can be calculated by,

$$\frac{\partial L}{\partial z'} = \frac{-u'}{z^d \sqrt{-f u'^2 - 2 u' z' +1}}\\ \frac{\partial L}{\partial z} = \frac{-u'^2 \partial_z f}{2 z^d \sqrt{-f u'^2 - 2 u' z' +1}}\\ \frac{\partial L}{\partial u'} = \frac{-f u' -z'}{z^d \sqrt{-f u'^2 - 2 u' z' +1}}\\ \frac{\partial L}{\partial u} = \frac{-u'^2 \partial_u f}{2 z^d \sqrt{-f u'^2 - 2 u' z' +1}}$$

where (I have written $f$ only to declutter)

$$f(z,u) = 1 - m(u) z^{d+1}, \qquad \frac{d}{dx} f = \left( \partial_z f \right) z' + \left( \partial_u f \right) u'$$

EOM for $z$: $\frac{d}{dx}\frac{\partial L}{\partial z'} - \frac{\partial L}{\partial z} = 0$

EOM for $u$: $\frac{d}{dx}\frac{\partial L}{\partial u'} - \frac{\partial L}{\partial u} = 0$

The boundary conditions (BC) are,

$z(x_0) = \epsilon, z'(x_s) = 0, u(x_0) = t, u'(x_s) = 0$

where the domain is $[x_0, x_s]$, $\epsilon$ is some arbitrary small number, and $t$ represents time for which I can choose some value.

I have not shown all the details, just the general picture of what I'm using. The details are in the code below.

d = 4;
x0 = 10^-5;
xs = 1;
c = 10^-2;
m0 = 1;
m[x_] := ((d + 1)/(c u[x] + (d + 1) m0^(-1/(d + 1))))^(d + 1);
f[x_] := 1 - m[x] z[x]^(d + 1);
fu[x_] := c z[x]^(d + 1) m[x]^((d + 2)/(d + 1));(*partial u*)
fz[x_] := -(d + 1) m[x] z[x]^d;(*partial z*)
fp[x_] := fz[x] z'[x] + fu[x] u'[x];(*x derivative of f*)
EOMz = -2 z[x]^d u''[x] - 2 z[x]^d u'[x]^2 z''[x] + (2 d f[x]^2 z[x]^(d - 1) - f[x] fz[x] z[x]^d) u'[x]^4 + (6 d f[x] z[x]^(d - 1) - 2 fz[x] z[x]^d) u'[x]^3 z'[x] + 4 d z[x]^(d - 1) u'[x]^2 z'[x]^2 - fp[x] z[x]^d u'[x]^3 + (fz[x] z[x]^d - 4 d f[x] z[x]^(d - 1)) u'[x]^2 - 6 d z[x]^(d - 1) u'[x] z'[x] + 2 d z[x]^(d - 1);
EOMu = -2 f[x] z[x]^d u''[x] + (2 z[x]^d u'[x] z'[x] - 2 z[x]^d) z''[x] - 2 d f[x]^2 z[x]^(d - 1) u'[x]^3 z'[x] - 4 d z[x]^(d - 1) u'[x] z'[x]^3 - 6 d f[x] z[x]^(d - 1) u'[x]^2 z'[x]^2 + fp[x] z[x]^d u'[x]^2 z'[x] + 2 d f[x] z[x]^(d - 1) u'[x] z'[x] + 2 d z[x]^(d - 1) z'[x]^2 - fp[x] z[x]^d u'[x];
s = NDSolveValue[Rationalize[{EOMz == 0, EOMu == 0, z[x0] == 10^-5, z'[xs] == 10^-16, u[x0] == 1, u'[xs] == 10^-16}, 0], {z, u}, {x, x0, xs}, Method -> {"Shooting", "StartingInitialConditions" -> {z[x0] == 10^-5, z'[x0] == 10^6, u[x0] == 1, u'[x0] == 10^-3}}, WorkingPrecision -> 20]
NDSolveValue::ndsz: At x == 0.99999810922483170542600892434185908076`20., step size is effectively zero; singularity or stiff system suspected.

Some comments about the code: I have chosen $10^{-16}$ for the vanishing derivatives at $x_s$ which is effectively zero, I chose $t = 1$. I employed shooting method for which I know that $z(x)$ should rapidly have a large value at the start so I guessed $z'(x_0) = 10^6$, for $u(x)$ I'm not sure of its initial behavior so I just guessed $u'(x_0) = 10^{-3}$. I have arbitrarily chosen some values for the initial mass $m_0$, the constant $c$, and the domain $[x_0, x_s]$.

The equations I have are coupled nonlinear second-order differential equations. As you see, Mathematica complains that the system of equations is stiff, I guess this can be alleviated by rescaling. However, I can first try to change the values of $m_0$ and $c$ to see if that will bring any changes, alas, none of it is working. Any advice on how to tackle this kind of problem?

*I don't suppose that Numerical GR techniques are needed in this case, however, I may be wrong.

Answer 1

We can solve this problem with using the Euler wavelets collocation method. First, we compute equations directly from L with EulerEquations as follows

Needs["VariationalMethods`"]
m = ((d + 1)/(c u[x] + (d + 1) m0^(-1/(d + 1))))^(d + 1); f = 
 1 - m z[x]^(d + 1); L = 
 Sqrt[-f u'[x]^2 - 2 u'[x] z'[x] + 1]/z[x]^d; eq1 = 
 EulerEquations[L, u[x], x];
eq2 = EulerEquations[L, z[x], x]; s = 
 Solve[{eq1, eq2}, {u''[x], z''[x]}][[1]] // Simplify;
eq01 = u''[x] == s[[1, 2]];
eq02 = z''[x] == s[[2, 2]];

Second, we define functions in wavelets basis and compute system of algebraic equations

UE[m_, t_] := EulerE[m, t];
psi[k_, n_, m_, t_] := 
  Piecewise[{{2^(k/2) UE[m, 2^k t - 2 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}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; zcol = 
 xcol = Table[(xl[[l - 1]] + xl[[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]; 
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
M = nn;
A = Array[a, {M}]; B = Array[b, M];
z2[x_] := A . Psi[x]; z1[x_] := A . int1[x] + a0; 
z0[x_] := A . int2[x] + a0 x + a1; u2[x_] := B . Psi[x]; 
u1[x_] := B . int1[x] + b0; u0[x_] := B . int2[x] + b0 x + b1;
var = Join[A, 
  B, {a0, a1, b0, b1}]; eqs = {eq01, eq02} /. {u''[x] -> u2[x], 
   u'[x] -> u1[x], u[x] -> u0[x], z''[x] -> z2[x], z'[x] -> z1[x], 
   z[x] -> z0[x]};
eq = Table[eqs /. {d -> 4, c -> 1/100, m0 -> 1}, {x, xcol}] // 
  Flatten; bc = {z0[x0] == x0, z1[xs] == 0, u0[x0] == 1, 
   u1[xs] == 0} /. {x0 -> 10^-5, xs -> 1};

Finally, we compute solution and plot u, z, L

sol = FindRoot[Join[eq, bc], 
  Table[{var[[i]], 9/10}, {i, Length[var]}], MaxIterations -> 1000];
{Plot[Evaluate[u0[x] /. sol], {x, 10^-5, 1}, AxesLabel -> {"x", "u"}, 
  PlotRange -> All], 
 Plot[Evaluate[z0[x] /. sol], {x, 10^-5, 1}, AxesLabel -> {"x", "z"}, 
  PlotRange -> All],LogPlot[L /. {u'[x] -> u1[x], u[x] -> u0[x], z'[x] -> z1[x], 
 z[x] -> z0[x]} /. {d -> 4, c -> 1/100, m0 -> 1} /. sol, {x, 

10^-5, 1}, PlotRange -> All, Frame -> True, PlotLabel -> "L"]}

Action

S = 
 NIntegrate[
  L /. {u'[x] -> u1[x], u[x] -> u0[x], z'[x] -> z1[x], 
      z[x] -> z0[x]} /. {d -> 4, c -> 1/100, m0 -> 1} /. sol, {x, 
   10^-5, 1}]
(* Out[]= 3.2155610^14 )

We can try to solve this problem with Haar wavelets as well as follows

Clear["Global`*"]
Needs["VariationalMethods`"]
m = ((d + 1)/(c u[x] + (d + 1) m0^(-1/(d + 1))))^(d + 1); f = 
 1 - m z[x]^(d + 1); L = 
 Sqrt[-f u'[x]^2 - 2 u'[x] z'[x] + 1]/z[x]^d; eq1 = 
 EulerEquations[L, u[x], x];
eq2 = EulerEquations[L, z[x], x]; s = 
 Solve[{eq1, eq2}, {u''[x], z''[x]}][[1]] // Simplify;
eq01 = u''[x] - s[[1, 2]] == 0;
eq02 = z''[x] - s[[2, 2]] == 0;
J = 5; M = 2^J; dx = 1/(2M); xl = Table[ldx, {l, 0, 2M}]; 
  xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, 2M + 1}]; 
  h1[x_] := Piecewise[{{1, 0 <= x <= 1}, {0, True}}]; 
  p1[x_, n_] := (1/n!)x^n; h[x_, k_, m_] := 
   Piecewise[{{1, Inequality[k/m, LessEqual, x, Less, 
       (1 + 2k)/(2m)]}, {-1, Inequality[(1 + 2k)/(2m), LessEqual, 
       x, Less, (1 + k)/m]}}, 0]
p[x_, k_, m_, n_] := Piecewise[{{0, x < k/m}, 
    {(-(k/m) + x)^n/n!, Inequality[k/m, LessEqual, x, Less, 
      (1 + 2k)/(2m)]}, 
    {((-(k/m) + x)^n - 2(-((1 + 2k)/(2m)) + x)^n)/n!, 
     (1 + 2k)/(2m) <= x <= (1 + k)/m}, 
    {((-(k/m) + x)^n + (-((1 + k)/m) + x)^n - 
       2(-((1 + 2k)/(2*m)) + x)^n)/n!, x > (1 + k)/m}}, 0]
var1 = Flatten[
  Table[a[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]]; var2 = 
 Flatten[Table[b[i, j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}]];
z2[x_] := 
 Sum[a[i, j] h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  a0 h1[x]; 
z1[x_] := 
 Sum[a[i, j] p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  a0 p1[x, 1] + a1; 
z0[x_] := 
 Sum[a[i, j] p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  a0 p1[x, 2] + a1 x + a2; 
u2[x_] := 
 Sum[b[i, j] h[x, i, 2^j], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  b0 h1[x]; 
u1[x_] := 
 Sum[b[i, j] p[x, i, 2^j, 1], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  b0 p1[x, 1] + b1; 
u0[x_] := 
 Sum[b[i, j] p[x, i, 2^j, 2], {j, 0, J, 1}, {i, 0, 2^j - 1, 1}] + 
  b0 p1[x, 2] + b1 x + b2;
eqs = {eq01, eq02} /. {u''[x] -> u2[x], u'[x] -> u1[x], u[x] -> u0[x],
     z''[x] -> z2[x], z'[x] -> z1[x], z[x] -> z0[x]};
var = Join[var1, var2, {a0, a1, a2, b0, b1, b2}];
eq = Table[eqs /. {d -> 4, c -> 1/100, m0 -> 1}, {x, xcol}] // 
  Flatten; bc = {z0[x0] == x0, z1[xs] == 0, u0[x0] == 1, 
   u1[xs] == 0} /. {x0 -> 10^-5, xs -> 1};
sol = FindRoot[Join[eq, bc], 
   Table[{var[[i]], 10^-1}, {i, Length[var]}], MaxIterations -> 1000, 
   Method -> {"Newton", "StepControl" -> "TrustRegion"}];
{Plot[Evaluate[u0[x] /. sol], {x, 10^-5, 1}, AxesLabel -> {"x", "u"}, 
  PlotRange -> All, PlotPoints -> 200], 
 Plot[Evaluate[z0[x] /. sol], {x, 10^-5, 1}, AxesLabel -> {"x", "z"}, 
  PlotRange -> All, PlotPoints -> 200], 
 LogPlot[L /. {u'[x] -> u1[x], u[x] -> u0[x], z'[x] -> z1[x], 
      z[x] -> z0[x]} /. {d -> 4, c -> 1/100, m0 -> 1} /. sol, {x, 
   10^-5, 1}, PlotRange -> All, Frame -> True, PlotLabel -> "L"]}

Action

S = 
 NIntegrate[
  L /. {u'[x] -> u1[x], u[x] -> u0[x], z'[x] -> z1[x], 
      z[x] -> z0[x]} /. {d -> 4, c -> 1/100, m0 -> 1} /. sol, {x, 
   10^-5, 1}, PrecisionGoal -> 4]
Out[]= 3.96352

Note, that action with Haar wavelets looks more optimal than that with Euler wavelets. This is probably a right solution. If we increased number of colocations points up to 128 (J=6), then the action increases to 5.60918.