I am trying to solve the coupled differential equation numerically with Mathematica. But the range of values are large so mathematica cannot give correct plot. Here is the code:

a = 4.75388*10^26;
b = 5.424*10^-3;
d = 4.75388*10^20;
{X, Y} = {x, y} /. 
z] == -((a/z) (x[z] - b*z^(3/2) E^(-z)) (BesselK[1, z]/BesselK[2, z])), 
y'[z] == ((d/z) (x[z] - b *z^(3/2) E^(-z)) (BesselK[1, z]/
BesselK[2, z]) - (a *z/4) (BesselK[1, z]) y[z]), 
x[0.1] == 1.552*10^-4, y[0.1] == 10^(-9)}, {x, y}, {z,0.1,100}] // 
FullSimplify // First
LogLogPlot[{X[z], Y[z]}, {z, 0.1, 100}, PlotRange -> All]

I am getting an Interpolating function whose domain is restricted to only {0.1,0.1}. I think this is the main reason of all the errors

  • 1
    $\begingroup$ This may be an obvious question: What about normalization of solution (to between 0 and 1) and then plotting it? Does that help? $\endgroup$ – dearN Sep 17 '19 at 11:48
  • $\begingroup$ When I run your code the first error I get says NDSolve::ndlim: Range specification z is not of the form {x, xend} or {x, xmin, xmax}. $\endgroup$ – N.J.Evans Sep 17 '19 at 12:23
  • $\begingroup$ Are we sure the issue is the size of the numbers and not the 5 errors that are given on evaluation? $\endgroup$ – ktm Sep 17 '19 at 12:57
  • $\begingroup$ When I run the code, NDSolve fails to take the first step (integration step); error NDSolve::ndsz. If it did succeed, the plot would fail because X is not a function, but a replacement Rule; hence X[z], while syntactically allowable, is semantically meaningless. If NDSolve can be fixed, try x[z] /. X as shown in the documentation. $\endgroup$ – Michael E2 Sep 18 '19 at 11:28

Geneal approach to problems like this is to attempt to regularize your DE. Naive approach in this case would be to replace $x$ with $x1=x/a$ and $y$ with $y1=y/d$:

a = 4.75388*10^26;
b = 5.424*10^-3;
d = 4.75388*10^20;

{X, Y} = {x1[#]*a &, y1[#]*d &} /. NDSolve[
     x1'[z] == -((BesselK[1, z] (-b E^-z z^(3/2) + a x1[z]))/(
       z BesselK[2, z]))
     , y1'[z] == (BesselK[1, z] (-b E^-z z^(3/2) + a x1[z]))/(
       z BesselK[2, z]) - (a z BesselK[1, z] y1[z])/4 
     , x1[0.1] == 1.552*10^-4/a
     , y1[0.1] == 10^(-9)/d
    , {x1, y1}
    , {z, 0.1, 100}
    ] // First
LogLogPlot[{X[z], Y[z]}, {z, 0.1, 100}, PlotRange -> All]


Resulting $y(z)$ does satisfy boundary conditions:

In[]:= Y[0.1]

Out[]= 1.*10^-9

But it might be not so obvious from picture above.

Plot[Y[z], {z, 0.1, 0.1 + 10^-10}, PlotRange -> All]

There is also an interval in which $x(z)$ becomes negative that is displayed somewhat misleading on the LogLogPlot

| improve this answer | |
  • 1
    $\begingroup$ it is a very approach but why not the y[z] initially starting from near the 10^-9 point since it is one of the given boundary conditions to solve this differential equation $\endgroup$ – user105697 Sep 18 '19 at 10:04
  • 1
    $\begingroup$ @user105697 edited reply to your concern into my answer. $\endgroup$ – Markhaim Sep 18 '19 at 10:36
  • $\begingroup$ Thank you for your effort. I think i have to revisit my input equations carefully. $\endgroup$ – user105697 Sep 18 '19 at 11:29
  • $\begingroup$ BTW, a more straightforward way to get the functions X, Y is to replace {x1[#]*a &, y1[#]*d &} /. NDSolve[args...] // First with NDSolveValue[args...]. $\endgroup$ – Michael E2 Sep 18 '19 at 11:48

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.