3
$\begingroup$

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.

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
3
$\begingroup$ 
$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

enter image description here

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

enter image description here

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

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks, did not notice that the issue was from numerical precision. $\endgroup$
    – Ali Hashmi
    Feb 25 at 1:29
3
$\begingroup$

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

Mathematica graphics

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]

Mathematica graphics

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]

Mathematica graphics

May be a bug.

Improve this answer
$\endgroup$
9
  • $\begingroup$ 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? $\endgroup$
    – Ali Hashmi
    Feb 25 at 1:03
  • $\begingroup$ @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 $\endgroup$
    – Nasser
    Feb 25 at 1:08
  • $\begingroup$ check this imgur.com/HE9i226 $\endgroup$
    – Ali Hashmi
    Feb 25 at 1:17
  • $\begingroup$ 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 $\endgroup$
    – Ali Hashmi
    Feb 25 at 1:21
  • $\begingroup$ 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. $\endgroup$
    – Ali Hashmi
    Feb 25 at 1:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.