Get a numerical solution to a nonlinear ODE?

Questions
, ,
1

I am trying to solve a nonlinear ODE BY applying a NDsolve and using StiffnessSwitching method, but when I try to find the root of my equation it gives me a error message. The same code is working well in Mathematica version 9, but in M.version 11.3 that I just upgraded is not working I do not know why this happened.
Would anyone help me please?

Here is my code

Z=800;
g= 0.023800000000000000000;
k2= 0.000194519;
 R= 1.5472;
ytest0= -13.911917694213733`;
? = $MachineEpsilon ;
y1[ytest_?NumericQ] := 
    NDSolve[
      {y''[r] + 2 y'[r]/r == ?2 Sinh[y[r]] , y[1] == ytest, 
       y'[?] == 0}, y, {r, ?, 1}, 
       Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
 y2[ytest_?NumericQ] := 
   NDSolve[
     {y''[r] + 2 y'[r]/r == ?2 Sinh[y[r]], 
      y[1] == ytest, y'[R] == 0}, y, {r, 1, R}, 
      Method -> {"StiffnessSwitching", "NonstiffTest" -> False}];
y1Try[ytest_?NumericQ] := y'[1] /. y1[ytest];
y2Try[ytest_?NumericQ] := y'[1] /. y2[ytest];
f = ytest /. FindRoot[y1Try[ytest] - y2Try[ytest]==-Zg, {ytest, ytest0}]
2

Hi,

I notice there is a typo in the code:

\[Kappa]2 = 0.000194519 (*it should be kappa instead of k*)

Also I highlighted ParametricNDSolveValue (doc) to simplify some parts of your code. The notebook with code is attached.

1685Sketch.png
1685Sketch.png1270×678 229 KB

ndsolve.nb (6.94 KB)

3

Why did you use (ParametricNDSolveValue)? Also, I did not see that you used ytest0.

4

The find root is not what I post it, when I write it as (see below) is not working??

f = ytest /. FindRoot[y1Try[ytest] - y2Try[ytest]==-Z/g, {ytest, ytest0}]
5

Hi,

ParametricNDSolveValue is basically the same as your y1 and y2, which are the functional of the solution of the two ODE’s. Looking at your code, I assume you want to know the 1st order derivative of each functional w.r.t. ytest parameter at location x=1. ParametricNDSolveValue kind of simplifies this process for you. Please correct me if I interpret your question incorrectly.

The other benefit of using ParametricNDSolveValue is that this function only processes the input equation once ( aka. conversion to canonical form) even the resulting functional varies as the parameters get updated.

6

You are right! Also, thank you for the great explanation about ParametricNDSolveValue. The question want to ask is why find root is not what I post it, when I write it as (see below) is not working??

f = ytest /. FindRoot[y1Try[ytest] - y2Try[ytest]==(-Zg), {ytest,
ytest0}] I did not see the discontinuity (Zg) in y’[r] at x == 1

7

I checked the equation to see where the zero might be through DiscretePlot from -5 to 5 with step=0.5

ScreenShot2018-09-17at9.11.03PM.png
ScreenShot2018-09-17at9.11.03PM.png517×284 9.59 KB

Also I notice y1Try[ytest] stops working when ytest<-11.5. Not sure about the mathematical reason behind this non-linear equation that causing the solver’s failure at this domain.

8

Hi

NDSolveValue[{y''[r] + 2 y'[r]/r == \[Kappa]2 Sinh[y[r]], y[1] == ytest0,
   y'[\[Epsilon]] == 0}, y, {r, \[Epsilon], 1}, 
 Method -> {"Shooting", 
   "StartingInitialConditions" -> {y[\[Epsilon]] == -11}}]

I found the shooting method extend the range of ytest0 to cover your input. -11 is found by using the StiffnessSwitching method first.