# NDSolve fails first order ODEs, but works when transformed to second order [duplicate]

The short version is that I am trying to solve a relatively simple system of coupled first order ODEs with NDSolve in Mathematica 12.1.1.0 on Windows. This system is perfectly linear in its unknowns (but depends non-linearly on the independent variable). Oddly, when I try to directly solve my boundary value problem, Mathematica fails (The error is NDSolve::nlnum1: The function value {...} is not a list of numbers with dimensions ... when the arguments are {...}.. However, when I eliminate one of the variables by hand (reducing the number of equations but increasing the differential order), it works. I'm puzzled because the system is, once again, linear, and this shouldn't happen!

In some more detail, I am trying to numerically solve a system of four coupled first order ODEs to obtain $$r(R)$$, $$g(R)$$, $$\Gamma(R)$$ and $$T_R(R)$$ (the subscript is not a derivative). This is a boundary value problem with bc $$r(A)=g(A)=g'(A)=A$$ and $$T_R(1)=A$$. Here, the domain goes from $$A$$ to $$1$$, where I expect the problem to be well-behaved for $$A\rightarrow0$$. The equation are linear in $$r$$,$$g$$,$$\Gamma$$,$$T_R$$ but depend non-linearly on $$R$$:

odes = {JgS^(2/3)
R^2 (3 JgS \[Chi] + 2 JgS (3 Rs \[Lambda] + \[Chi]) +
6 JgS TR[R]) + (12 JgS \[Lambda] +
JgS^(5/3) R^2 \[Chi]) \[CapitalGamma][R] ==
2 (6 JgS \[Lambda] +
JgS^(4/3) R \[Chi] (2 g[R] + R Derivative[1][g][R])),
JgS^(1/3) (-4 (3 JgS \[Mu] + JgS \[Chi]) g[R] + 12 JgS \[Mu] r[R] +
R (JgS^(
1/3) (3 JgS \[Chi] +
2 JgS (3 Rs \[Lambda] + 3 \[Mu] + \[Chi])) +
6 JgS^(4/3)
TR[R] + (-6 JgS \[Mu] +
JgS \[Chi]) (JgS^(1/3) \[CapitalGamma][R] -
2 Derivative[1][g][R]) -
12 JgS \[Mu] Derivative[1][r][R])) +
6 JgS \[Lambda] Derivative[1][\[CapitalGamma]][R] == 0,
4 \[Mu] g[R] +
R (2 JgS^(1/3) \[Mu] \[CapitalGamma][R] +
4 \[Mu] (-Derivative[1][g][R] + Derivative[1][r][R]) +
JgS^(1/3) (-2 \[Mu] + R Derivative[1][TR][R])) ==
4 \[Mu] r[R],
4 \[Kappa] g[R] + JgS^(1/3) R (\[Kappa] + 2 \[Mu] + 2 TR[R]) ==
4 \[Kappa] r[R] +
R (\[Kappa] + 2 \[Mu]) (JgS^(1/3) \[CapitalGamma][R] -
2 Derivative[1][g][R] + 2 Derivative[1][r][R])};


Here is an attempt to solve this with NDSolve, for some parameters:

A = 10^-2;
subParam = {\[Mu] -> 1, \[Lambda] -> 5, \[Chi] -> 0.1,
JgS -> 5, \[Kappa] -> 1, Rs -> 0.3};
sol = NDSolve[{odes, r[A] == 0, TR[1] == 0,  g'[A] == A,
g[A] == A} /. subParam, {r[R], g[R],
TR[R], \[CapitalGamma][R]}, {R, A, 1}];


This produces an error message NDSolve::nlnum1: The function value {0.,3450.87,-(1/100)+(g^\[Prime])[1/100],-0.01} is not a list of numbers with dimensions {4} when the arguments are {0.,0.,0.,0.,-586608.,-583596.,3450.87,-11008.}.

Now, let us try another approach: We use the first equation to eliminate $$\Gamma$$ by hand:

First@Solve[odes[[1]], \[CapitalGamma][R]];
sub\[CapitalGamma] =
Head[%[[1, 1]]] -> Function[R, Evaluate[%[[1, 2]]]]
temp = odes[[2 ;; -1]] /. sub\[CapitalGamma];


The variable temp now contains the set of reduced ODEs. The number of equations and unknowns have been reduced by one, but the differential order of $$g$$ has increased. So all should work just the same. Now let us numerically solve this and plot the solution:

   A = 10^-3;
subParam = {\[Mu] -> 1, \[Lambda] -> 5, \[Chi] -> 0.1,
JgS -> 5, \[Kappa] -> 1, Rs -> 0.3};
sol = NDSolve[{temp, r[A] == 0, TR[1] == 0,  g'[A] == A,
g[A] == A} /. subParam, {r[R], g[R], TR[R], g'[R](*
needed for \[CapitalGamma][R]*)}, {R, A, 1}];
Plot[{r[R], g[R],  \[CapitalGamma][R], TR[R] } /.
sub\[CapitalGamma] /. Jg -> JgS /. subParam /. sol //
Evaluate, {R, A, 1},
PlotLegends -> {"r", "g", "\!$$\*SubscriptBox[\(T$$, $$R$$]\)",
"\[CapitalGamma]" }, PlotTheme -> "Scientific",
PlotRange -> {All, {-0.3, 1.2}}]


The reduced equation can be solved as expected! How can this be? Between the full and reduced system lies just a linear transformation.

The motivation for my question is that here we're dealing with a linearised system; the fully non-linear one is implicit, and I can't symbolically solve for $$\Gamma(R)$$ as I have done here.

• Have you seen this question and the answere there? System of ODEs causes NDSolve::nlnum1 warning Commented Apr 29 at 13:26
• Thanks. The reply in that question hints that maybe my boundary condition $g'(A)=A$ may be the reason, but I have no mixed BC like in that question. Besides, my BC corresponds to $\Gamma=1$ in the limit $A\rightarrow0$ so I'm not sure where I'm being unreasonable. Commented Apr 30 at 7:32
• In general, the derivative order of the boundary conditions should be less than the derivative order of the differential equations. Isolate g'[A] by Equal @@ Simplify[Solve[odes, {r'[R], g'[R], TR'[R], \[CapitalGamma]'[R]}][[1, 2]]] and then use Simplify[% /. R -> A /. {r[A] -> 0, TR[1] -> 0, g'[A] -> A, g[A] -> A}] to provide a new boundary condition to replace the one involving g'[A] . Commented May 1 at 3:28
• "but I have no mixed BC like in that question." Please see the simplified example in that answer of mine, there's nothing to do with mixed boundary condition. The only related part is, as pointed out by other comments and answer, the differential order of the boundary condition $g'(A)=A$ is too high. Commented May 1 at 4:15

To illustrate my comment, obtain a formula for g'[R] from odes

Equal @@ Simplify[Solve[odes, {r'[R], g'[R], TR'[R], \[CapitalGamma]'[R]}][[1, 2]]];


apply the R = A boundary conditions to it

Simplify[% /. R -> A /. {r[A] -> 0, TR[1] -> 0, g'[A] -> A, g[A] -> A}]


and use the result in place of the g'[R] boundary condition.

sol = NDSolveValue[{odes, r[A] == 0, TR[1] == 0, %, g[A] == A} /.
subParam, {r[R], g[R], TR[R], \[CapitalGamma][R]}, {R, A, 1}] // Flatten;

Plot[%, {R, A, 1}, AxesLabel -> {R, "r,g,TR,\[CapitalGamma]"},
LabelStyle -> {12, Bold, Black}


• Thank you, this makes sense! But is there any way this solution could be sped up? It takes very long compared to my proposition to eliminate one of the variables analytically. I have tried your approach on a more complex system (non-linear version) and aborted the solution after ~5min of waiting... Commented May 1 at 10:19
• @AlexanderErlich, in general the speed of NDSolve depends on the details of the ODEs. But here are two guesses. \[CapitalGamma][R] varies rapidly for small R, causing the NDSolve step size in R to be very small. So, eliminating it from the ODEs may increase the step size substantially. Another possibility is that eliminating \[CapitalGamma][R] allows the boundary value computation to converge more rapidly. Try eliminating [CapitalGamma][R]` as you have done but with my modified boundary conditions and see what happens. Commented May 1 at 15:01