Message Boards Message Boards

Get a numerical solution to a nonlinear ODE?

Posted 6 years ago

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}]
POSTED BY: kolod al
7 Replies

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.

pndsolvev

Attachments:
POSTED BY: Shenghui Yang
Posted 6 years ago

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

POSTED BY: kolod al
Posted 6 years 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}]
POSTED BY: kolod al

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.

POSTED BY: Shenghui Yang
Posted 6 years 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 BY: kolod al

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

plot

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.

POSTED BY: Shenghui Yang

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.

POSTED BY: Shenghui Yang
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract