# Why does introducing value of variable for NDSolve via function not work for this problem?

I am very confused and even restarted Kernel several times (Mathematica 13.1 Windows 10). However, I fail to understand why one evaluation of NDSolve works while a similar implementation in which I introduce the value of a scalar variable LK into the equation via function definition does not work.

This works

dP = 17.48;
c0 = -4.0;
LK = -39.9968;

First@NDSolve[{c'[s] - (-2 d[s] Sqrt[1 - f[s] c[s]^2])/f[s] == 0,
f'[s] - 4 Sqrt[1 - f[s] c[s]^2] == 0,
d'[s] - (-(2 c[s]^2 (d[s] - c0) + c[s] (c0^2 - d[s]^2) + LK*c[s] +
dP + 4 d[s] (1 - f[s] c[s]^2)/f[s]))/(Sqrt[
1 - f[s] c[s]^2]) == 0, f[0.0005] == 0.0001,
c[0.0005] == -0.50, d[0.0005] == 0.0001}, {c, f, d}, {s, 0.0005, 0.5}]


But this does not

dP = 17.48;
c0 = -4.0;

psol[LK_?NumericQ] :=
First@NDSolve[{c'[s] - (-2 d[s] Sqrt[1 - f[s] c[s]^2])/f[s] == 0,
f'[s] - 4 Sqrt[1 - f[s] c[s]^2] == 0,
d'[s] - (-(2 c[s]^2 (d[s] - c0) + c[s] (c0^2 - d[s]^2) + LK*c[s] +
dP + 4 d[s] (1 - f[s] c[s]^2)/f[s]))/(Sqrt[
1 - f[s] c[s]^2]) == 0, f[0.0005] == 0.0001,
c[0.0005] == -0.50, d[0.0005] == 0.0001}, {c, f, d}, {s, 0.0005, 0.5}]

psol[-33.9968]


Does anyone experience the same issue? What am I possibly doing wrong here.

\$Version

(* "13.2.1 for Mac OS X ARM (64-bit) (January 27, 2023)" *)


Your second method evaluates very much more slowly, stops evaluating at about s = 0.433, and gives erroneous results.

Clear["Global*"]

dP = 17.48;
c0 = -4.0;

psol[LK_?NumericQ] :=
First@NDSolve[{c'[s] - (-2 d[s] Sqrt[1 - f[s] c[s]^2])/f[s] == 0,
f'[s] - 4 Sqrt[1 - f[s] c[s]^2] == 0,
d'[s] - (-(2 c[s]^2 (d[s] - c0) + c[s] (c0^2 - d[s]^2) + LK*c[s] + dP +
4 d[s] (1 - f[s] c[s]^2)/f[s]))/(Sqrt[1 - f[s] c[s]^2]) == 0,
f[0.0005] == 0.0001, c[0.0005] == -0.50, d[0.0005] == 0.0001}, {c, f,
d}, {s, 0.0005, 0.5}]

psol[-33.9968] // AbsoluteTiming


Assuming that this reflects a precision issue, increase the precision handling by using exact values for all constants.

Clear["Global*"]

dP = 1748/100;
c0 = -4;

psol[LKv_?NumericQ] :=
Module[{LK = Rationalize[LKv, 0]},
NDSolve[{c'[s] - (-2 d[s] Sqrt[1 - f[s] c[s]^2])/f[s] == 0,
f'[s] - 4 Sqrt[1 - f[s] c[s]^2] == 0,
d'[s] - (-(2 c[s]^2 (d[s] - c0) + c[s] (c0^2 - d[s]^2) + LK*c[s] + dP +
4 d[s] (1 - f[s] c[s]^2)/f[s]))/(Sqrt[1 - f[s] c[s]^2]) == 0,
f[5*^-4] == 10^-4, c[5*^-4] == -1/2, d[5*^-4] == 10^-4}, {c, f, d}, {s,
5*^-4, 1/2}] // First]

psol[-39.9968] // AbsoluteTiming


Plot[Evaluate[{c[s], f[s], d[s]} /. psol[-39.9968]],
{s, 0.0005, 0.5},
PlotLabels -> {c[s], f[s], d[s]}]


• Thanks, did not notice that the issue was from numerical precision. Feb 25 at 1:29

it looks like a naming/scoping/buffering issue if you ask me. I can't answer why, but a simple workaround is the following. This is something that should be done anyway. Which is to define the ode and the IC as separate variables.

This makes the code easier to read also and maintain, and resolves this issue.

dP = 17.48;
c0 = -4.0;
LK = -39.9968;
ode = {c'[s] - (-2 d[s] Sqrt[1 - f[s] c[s]^2])/f[s] == 0,
f'[s] - 4 Sqrt[1 - f[s] c[s]^2] == 0,
d'[s] - (-(2 c[s]^2 (d[s] - c0) + c[s] (c0^2 - d[s]^2) + LK*c[s] +
dP + 4 d[s] (1 - f[s] c[s]^2)/f[s]))/(Sqrt[
1 - f[s] c[s]^2]) == 0};
ic = {f[0.0005] == 0.0001, c[0.0005] == -0.50, d[0.0005] == 0.0001};

NDSolve[{ode, ic}, {c, f, d}, {s, 0.0005, 0.5}]


Now the second one works with no problem

dP = 17.48;
c0 = -4.0;
psol[LK_?NumericQ] := First@NDSolve[{ode, ic}, {c, f, d}, {s, 0.0005, 0.5}]

psol[-33.9968]


Without doing the above, it still works, but takes much longer and gives

dP = 17.48;
c0 = -4.0;
psol[LK_?NumericQ] :=
First@NDSolve[{c'[s] - (-2 d[s] Sqrt[1 - f[s] c[s]^2])/f[s] == 0,
f'[s] - 4 Sqrt[1 - f[s] c[s]^2] == 0,
d'[s] - (-(2 c[s]^2 (d[s] - c0) + c[s] (c0^2 - d[s]^2) + LK*c[s] +
dP + 4 d[s] (1 - f[s] c[s]^2)/f[s]))/(Sqrt[
1 - f[s] c[s]^2]) == 0, f[0.0005] == 0.0001,
c[0.0005] == -0.50, d[0.0005] == 0.0001}, {c, f, d}, {s, 0.0005,
0.5}]

psol[-33.9968]


May be a bug.

• I think if you do Remove[Lk] before your code then it does not work... dP = 17.48; c0 = -4.0; psol[LK_?NumericQ] := First@NDSolve[{ode, ic}, {c, f, d}, {s, 0.0005, 0.5}] psol[-33.9968] ... can you confirm? Feb 25 at 1:03
• @AliHashmi it works just fine. Why would removing LK causes a problem? Here is screen shot fyi !Mathematica graphics This is using 13.1 on windows 10 Feb 25 at 1:08
• check this imgur.com/HE9i226 Feb 25 at 1:17
• I just did ClearAll and Remove then added your code but did not define LK like the first time you defined it. I just cant seem to run the code. If you look at the output, LK stays undefined. I dont know what is happening on my end Feb 25 at 1:21
• I think you need to redefine ode after clearing or removing LK. because the ode still retains the value of -33.66 even after you clear variable LK. Feb 25 at 1:25