# NDSolve problem

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

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

## 2 Answers

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


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.

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