ParametricNDSolveValue Doesn't Work

Question

I am plotting particle trajectory using ParametricNDSolveValue, but it gives error. Any help is highly appreciated.

Clear["Global`*"]
h=1-rs/r+(8 \[Pi] r0^3 \[Rho]0)/Sqrt[r^2+r0^2]-(8 \[Pi] r0^3 \[Rho]0  ArcSinh[r/r0])/r;
g=1-rs/r[\[Phi]]-(8 \[Pi] r0^3 \[Rho]0 ArcSinh[r[\[Phi]]/r0])/r[\[Phi]]+(8 \[Pi] r0^3 \    
[Rho]0)/Sqrt[r0^2+r[\[Phi]]^2];
G=Assuming[u[\[Phi]]>0,1/b^2-u[\[Phi]]^2(g)/.r[\[Phi]]->1/u[\[Phi]]];
\[Phi]tailcontrolsigmoid[\[Phi]min_,\[Phi]max_,\[Sigma]_]:=(1/(1+Exp[(\[Phi]min/\[Phi]max-1)/\ 
[Sigma]]))
rs=2;
r0=1;
\[Rho]0=0.01;
bmax=1.2;
db=0.2;
\[Phi]min=0.001;
\[Phi]max=3\[Pi];
PhiMax[b_]:=(1-b/(bmax)) \[Phi]max
rF=ParametricNDSolveValue[{u''[\[Phi]]==1/2 D[G,u[\[Phi]]],u'[0]==1/b,u[0]==0},u,{\[Phi],\ 
[Phi]min,\[Phi]max*\[Phi]tailcontrolsigmoid[\[Phi]min,\[Phi]max,0.2]},b,StartingStepSize- 
>1/1000,Method->"ExplicitEuler"];
Horizon=FindRoot[h,{r,rs}][[1,2]];
BH=PolarPlot[{r=Horizon},{\[Theta],0,2\[Pi]},PlotStyle->{Directive[Dashed,Thick,Red]}];
BH1=Graphics[{PointSize[0.08],Point[{0,0}]}];
A3=Legended[Show@Table[PolarPlot[Style[1/Max[0.001,rF[b][\[Phi]]],Blend[{Cyan,Magenta}, 
(b)/(bmax)]],{\[Phi],\[Phi]min,rF[b]["Coordinates"][[1,-1]]},PlotRange->{{-bmax,bmax},{- 
bmax,bmax}},Frame->True,ImageSize->500], 
{b,1,bmax,db}],Placed[BarLegend[{(Blend[{Cyan,Magenta},(#-1)/(bmax-1)]&), 
{1,bmax}},LegendLayout->"Column",LegendLabel->"b"],Right]];
B=Show[A3,BH,BH1]

There are several typos in your code. Please, explain your model first. — Alex Trounev
– Alex Trounev, Commented Nov 1 at 11:10

Alex Trounev · Accepted Answer · 2025-11-01 12:32:55Z

I simplified the model to make it more understandable.

Clear["Global`*"]
h = 1 - rs/r + (8 \[Pi] r0^3 \[Rho]0)/
    Sqrt[r^2 + r0^2] - (8 \[Pi] r0^3 \[Rho]0 ArcSinh[r/r0])/r;
g = 1 - rs/r[\[Phi]] - (8 \[Pi] r0^3 \[Rho]0 ArcSinh[r[\[Phi]]/r0])/
    r[\[Phi]] + (8 \[Pi] r0^3 \[Rho]0)/Sqrt[r0^2 + r[\[Phi]]^2];
G = Assuming[u[\[Phi]] > 0, 
   1/b^2 - u[\[Phi]]^2 (g) /. r[\[Phi]] -> 1/u[\[Phi]]];
\[Phi]tailcontrolsigmoid[\[Phi]min_, \[Phi]max_, \[Sigma]_] := (1/(1 +
      Exp[(\[Phi]min/\[Phi]max - 1)/\[Sigma]]))
rs = 2;
r0 = 1;
\[Rho]0 = 0.01;
bmax = 1.2;
db = 0.2;
\[Phi]min = 0.001;
\[Phi]max = 3 \[Pi];
PhiMax[b_] := (1 - b/(bmax)) \[Phi]max
rF = ParametricNDSolveValue[{u''[\[Phi]] == 
     1/2 (-((-2 - 0.25132741228718347` ArcSinh[1/u[\[Phi]]] + 
             0.25132741228718347`/((1 + 1/u[\[Phi]]^2)^(3/2)
               u[\[Phi]]^3) + 0.25132741228718347`/(
             Sqrt[1 + 1/u[\[Phi]]^2] u[\[Phi]])) u[\[Phi]]^2) - 
        2 u[\[Phi]] (1 + 0.25132741228718347`/Sqrt[
           1 + 1/u[\[Phi]]^2] - 2 u[\[Phi]] - 
           0.25132741228718347` ArcSinh[1/u[\[Phi]]] u[\[Phi]])), 
    u'[0.001] == 1/b, u[0.001] == .001}, u, {\[Phi], .001, 2.4}, b];
Horizon = FindRoot[h, {r, rs}][[1, 2]];

Visualization

pl = Table[
  PolarPlot[1/rF[b][\[Phi]], {\[Phi], \[Phi]min, 2.4}, 
   PlotRange -> Automatic, Frame -> True, AspectRatio -> Automatic, 
   PlotStyle -> Blend[{Cyan, Magenta}, (b - 1)/(bmax - 1)]], {b, 1, 
   bmax, db}];

BH = PolarPlot[{r = Horizon}, {\[Theta], 0, 2 \[Pi]}, 
   PlotStyle -> {Directive[Dashed, Thick, Red]}];
BH1 = Graphics[{PointSize[0.08], Point[{0, 0}]}];
A3 = Legended[Show@pl, 
   Placed[BarLegend[{(Blend[{Cyan, 
          Magenta}, (# - 1)/(bmax - 1)] &), {1, bmax}}, 
     LegendLayout -> "Column", LegendLabel -> "b"], Right]];
B = Show[A3, BH, BH1]

This is more advanced code for bmax=4

rF = ParametricNDSolveValue[{u''[\[Phi]] == 
     1/2 (-((-2 - 0.25132741228718347` ArcSinh[1/u[\[Phi]]] + 
             0.25132741228718347`/((1 + 1/u[\[Phi]]^2)^(3/2)
               u[\[Phi]]^3) + 0.25132741228718347`/(
             Sqrt[1 + 1/u[\[Phi]]^2] u[\[Phi]])) u[\[Phi]]^2) - 
        2 u[\[Phi]] (1 + 0.25132741228718347`/Sqrt[
           1 + 1/u[\[Phi]]^2] - 2 u[\[Phi]] - 
           0.25132741228718347` ArcSinh[1/u[\[Phi]]] u[\[Phi]])), 
    u'[0.001] == 1/b, u[0.001] == .001}, 
   u, {\[Phi], .001, phi}, {b, phi}];

Visualization

With[{bmax = 4}, m = Table[2.4 + (x - 1)/2.4, {x, 1, bmax, db}]; 
 bb = Table[x, {x, 1, bmax, db}]; 
 pl1 = Table[
   PolarPlot[1/rF[bb[[i]], m[[i]]][\[Phi]], {\[Phi], \[Phi]min, m[[i]]},
     PlotRange -> Automatic, Frame -> True, AspectRatio -> Automatic, 
    PlotStyle -> 
     Blend[{Cyan, Magenta}, (bb[[i]] - 1)/(bmax - 1)]], {i, Length@bb}];
 BH = PolarPlot[{r = Horizon}, {\[Theta], 0, 2 \[Pi]}, 
   PlotStyle -> {Directive[Dashed, Thick, Red]}];
 BH1 = Graphics[{PointSize[0.08], Point[{0, 0}]}];
 A3 = Legended[Show@pl1, 
   Placed[BarLegend[{(Blend[{Cyan, 
          Magenta}, (# - 1)/(bmax - 1)] &), {1, bmax}}, 
     LegendLayout -> "Column", LegendLabel -> "b"], Right]];
 B = Show[A3, BH, BH1]]

Stack Exchange Network

ParametricNDSolveValue Doesn't Work

1 Answer 1

Your Answer

Hot Network Questions

ParametricNDSolveValue Doesn't Work

1 Answer 1

Your Answer

Sign up or log in

Post as a guest

Related

Hot Network Questions