Message Boards Message Boards


Get a numerical solution to a nonlinear ODE?

Posted 1 year ago
7 Replies
1 Total Likes

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

g= 0.023800000000000000000;
k2= 0.000194519;
 R= 1.5472;
ytest0= -13.911917694213733`;
? = $MachineEpsilon ;
y1[ytest_?NumericQ] := 
      {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] := 
     {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}]
7 Replies


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.


Posted 1 year ago

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


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.

Posted 1 year ago

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

Posted 1 year ago

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

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


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.


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.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract