2
$\begingroup$

I am searching for a while now, but I don't seem to be able to find an Answer for my Problem - if I am just not able to search properly, I am really sorry.

I Simplified my Problem to the following example:

eqnA := {a''[x] + a[x] == b*Sin[a[x]]};
condA := {a[0] == a'[0] == 1};
stopCond := {WhenEvent[Evaluate[Re[a[x]] <= 0], xMax = x; "StopIntegration"]}
system := Join[eqnA, condA, stopCond];
sol = ParametricNDSolve[system, a, {x, 0, \[Infinity]}, {b}]
Plot[Evaluate@Table[a[b][x] /. sol, {b, .8, 1.2, .05}], {x, 0, xMax},PlotRange -> All]

Now, the Plot I want to create is a Parametric Plot of a over xMax. (Yes, I understand, that in this example it would be just a line - if it is any help: I am trying to solve the TOV equations and create the Mass over Radius diagram with central density as parameter)

I think, my understanding of WhenEvent and the uses of it might be flawed.

Any Tips are appreciated!

Edit: Since I think I did not formulate my problem clearly enough I try to explain a bit further:

With the post from george2079 I was able to get a bit further, but It still does not work.

f[b_] := (a[b] /. sol)["Domain"][[1, 2]]
g[b_] := a[b][f[b]] /. sol
ParametricPlot[{f[b], g[b]}, {b, 0, 2}, PlotRange -> All]

I would expect that to produce something close to a line along the x-axis, but instead I get an error message. Plotting f[b] and g[b] individually works just fine.

$\endgroup$
  • $\begingroup$ By "Parametric Plot of a over xMax", do you mean a plot of the point {xMax, a[b][xMax]} as b varies? Note that a[b][xMax] would approximately be zero. $\endgroup$ – Michael E2 May 17 '16 at 14:09
  • $\begingroup$ If I understand you correctly, yes. xMax[b] on the x-axis, on the y-axis a[b][xMax], and I am interested in the resulting line when varying b (in my simple example a straight line along the x-axis) $\endgroup$ – Zoldor May 17 '16 at 14:16
3
$\begingroup$

Your WhenEvent is just fine. If you look at an individual solution you can see where it stopped:

            (a[.8] /. sol)["Domain"]

{{0., 3.6351}}

The trouble is Plot extrapolates* past the solution domain, so you need to specify each plot domain individually:

Show[Table[
  Plot[ #[x] /. sol ,
    Evaluate[{x, Sequence @@ (First@#["Domain"])}],
    PlotStyle -> Hue[RandomReal[]]] &@(a[b] /. sol) ,
     {b, .8, 1.2, .05} ], PlotRange -> All]

enter image description here

Of course knowing whats going on you might just specify PlotRange->{0,Automatic} as well.

If you want to plot the stop point vs. b :

Plot[ (a[b] /. sol)["Domain"][[1, 2]] , {b, .8, 1.2}]

enter image description here

*and it really ought to throw a warning about it...

Edit: here is a parametric plot of the end value of each solution:

ParametricPlot[
  (Flatten@Last@Transpose[(a[b] /. sol)[{"Grid", "ValuesOnGrid"}]]),
       {b, .8, 1.2}, PlotRange -> {Automatic, {-1, 1}}, 
       AxesOrigin -> {Automatic, -1}, AspectRatio -> 1/GoldenRatio]

enter image description here

letting ParametricPlot autorange we can see the accuracy of the solver meeting your zero condition.

enter image description here

$\endgroup$
  • $\begingroup$ Thank you, this already looks really useful! Maybe I can solve my Problem playing with "Domain". But I think my question formulated a bit confusing (Sorry for that, english is not my native tongue). I tried to explain, that the Plot I really want to create would be something like: ParametricPlot[{xMax[b],a[b]},{b, .8, 1.2}] $\endgroup$ – Zoldor May 18 '16 at 5:41
  • $\begingroup$ responding to comment, a[b] at what x? a[b][xMax[b]] is zero due to your WhenEvent $\endgroup$ – george2079 May 18 '16 at 15:25
  • $\begingroup$ yes, a[b][xMax[b]] is zero in my example code, but in the Parametric Plot it should be a Line on the x-axis (since a is on the y-axis and because of my condition). Like I said, this is because i used a simplified ODE as example. Since this causes a lot of confusion: My full system are the coupled Tolman-Oppenheimer-Volkoff equations, with the WhenEvent triggered by the pressure reaching approx. 0 and the Plot should be the Mass at r=rMax over the max. of the Radius (which I get via the WhenEvent). But the problem is the same as in the example ODE that I used. $\endgroup$ – Zoldor May 19 '16 at 4:40
  • $\begingroup$ Thank you, this does exactly what I need! $\endgroup$ – Zoldor May 20 '16 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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