1
$\begingroup$

I am trying to solve numerically this differential equation

s := NDSolve[{y''[x] + ω[x]*y[x] - 1/(y[x])^3 == 0,  y[0] == 1, y'[0] == 0.3}, y, {x, 0, 10}]

where ω[x] is a complicated function, given by:

ω[x]=E^(4 (1 + 1/3 Log[Sinh[0.0411881 x]] - (
   Cos[π/18] Log[(
     Sqrt[Cos[π/18]^2 + Sin[π/18]^2] + Sqrt[
      Cos[π/18]^2 + Sin[π/18]^2 + 
       0.00001 Sinh[0.0411881 x]^2])/(
     Sqrt[Cos[π/18]^2 + Sin[π/18]^2] - Sqrt[
      Cos[π/18]^2 + Sin[π/18]^2 + 
       0.00001 Sinh[0.0411881 x]^2])])/(
   6 Sqrt[Cos[π/18]^2 + Sin[π/18]^2])))

Since this function has a divergence in 0, Mathematica can't solve the equation above. It gives the following errors: Power::infy: Infinite expression 1/0. encountered. NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0..

Can you suggest any other way I can find the numerical solution?

I tried to move the interval, starting from x==0.001, but the solution obtained doesn't appear to be the right one.

Thank you all

$\endgroup$
  • $\begingroup$ Is there any reason that you're not simplifying your expression for ω[x]? Note that Sqrt[Cos[π/18]^2 + Sin[π/18]^2] = 1. $\endgroup$ – Michael Seifert Aug 2 '17 at 14:49
2
$\begingroup$

It always pays to simplify your equations:

ω2[x_] = Simplify[ω[x]]

(* E^4 Sinh[0.0411881 x]^(4/3) (-((1 + Sqrt[1. + 0.00001 Sinh[0.0411881 x]^2])/(-1 + Sqrt[1. + 0.00001 Sinh[0.0411881 x]^2])))^(-(2/3) Cos[π/18]) *)

Since this function is not indeterminate as $x \to 0$, you can use it in your ODE without a problem:

s = NDSolve[{y''[x] + ω2[x]*y[x] - 1/(y[x])^3 == 0, y[0] == 1, y'[0] == 0.3}, y, {x, 0, 10}]
Plot[Evaluate[ReIm[y[x] /. s]], {x, 0, 10}]

enter image description here

As @zhk noted in their answer, the function $\omega(x)$ is complex as defined, and so the resulting y[x] is complex as well. In the graph above, the real part of y[x] is in blue and the imaginary part is in yellow. The imaginary part stays small for the range of $t$ in the above plot, but it starts to become comparable to the real part around $t = 30$ or so.

$\endgroup$
1
$\begingroup$

The proper way to write a set is ω[x_] = not ω[x] =

ω[x_] = 
 E^(4 (1 + 
      1/3 Log[Sinh[
         0.0411881 x]] - (Cos[π/
           18] Log[(Sqrt[Cos[π/18]^2 + Sin[π/18]^2] + 
             Sqrt[Cos[π/18]^2 + Sin[π/18]^2 + 
               0.00001 Sinh[0.0411881 x]^2])/(Sqrt[
              Cos[π/18]^2 + Sin[π/18]^2] - 
             Sqrt[Cos[π/18]^2 + Sin[π/18]^2 + 
               0.00001 Sinh[0.0411881 x]^2])])/(6 Sqrt[
          Cos[π/18]^2 + Sin[π/18]^2])))

Now, if you assign a random value to x to check ω,

ω[1]

-1.16889*10^-6 - 2.18208*10^-6 I

the output shows that ω is complex. So, we need to use Re@ω if interested in real values only.

s = NDSolve[{y''[x] + Re@ω[x]*y[x] - 1/(y[x])^3 == 0, 
       y[10^-4] == 1, y'[10^-4] == 0.3}, y, {x, 10^-4, 10}]

Plot[y[x] /. s, {x, 10^-4, 10}]

enter image description here

$\endgroup$

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.