Numerically solving an ODE whose right-hand side involves an integral
Being a mathematician, I resist fudging by cutting off the singularity by some small eps = 10^-12
. But if you're an engineer, you should be satisfied @Nasser's answer, which likely yields single-precision accuracy, which is pretty good after all.
A little analysis shows
$$\eqalign{
w''(x) &= -{1\over x^2} \int_0^x \sin w(t) \; dt + {\sin w(x) \over x} \cr
&= \int_0^x {\sin w(x) - \sin w(t) \over x^2} \; dt \cr
&= \int_0^x {(x-t) \over x^2}{\sin w(x) - \sin w(t) \over x-t} \; dt \cr
&\buildrel \rm{MVT} \over = \int_0^x {(x-t) \over x^2}\,\cos w(\xi) \, w'(\xi)\,(x-t) \; dt \cr
&= \int_0^x {\left(1-{t\over x}\right)^2}\,\cos w(\xi) \, w'(\xi) \; dt \cr
}$$
where $\xi = \xi_t$ is between $x$ and $t$ and depends on $t$ for a fixed $x$. Of the three factors in the integrand as $x \rightarrow 0$, the first lies between $0$ and $1$, the cosine approaches zero since $w(\xi) \rightarrow \pi/2$, and $w'(\xi)$ approaches $1$. Therefore $w''(x) \rightarrow 0$ as $x \rightarrow 0$.
Thus we should define our ODE so that $w''(0) = 0$. Piecewise
is the tool for that.
sol = NDSolve[{w''[x] == Piecewise[{{(Sin[w[x]] - w'[x])/x, x != 0}}],
w[0] == Pi/2, w'[0] == 1}, w, {x, 0, 4 Pi}];
Plot[w[x] /. sol, {x, 0, 4 Pi}, AxesOrigin -> {0, 0}]
It differs from Nasser's answer by less than 3 * 10^-8
over the interval.
There is a singularity at x=0
.
First, differentiate and use Leibniz rule, this gives the non-linear ODE
$$ x w''(x)+ w'(x) = \sin(w(x)) $$
But NDSolve
can't solve this due to singularity at $x=0$
eq = x w''[x] + w'[x] == Sin[w[x]];
NDSolve[{eq, w[0] == Pi/2, w'[0] == 1}, w, {x, 0, 1}]
NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.
So using trick to bypass $x=0$ (maybe a better trick exists, but the idea is to avoid zero)
sol = With[{eps = 10^-12},
NDSolve[{w''[x] ==
If[x < eps, -w'[x]/eps + Sin[w[x]]/eps, -w'[x]/x + Sin[w[x]]/x],
w[0] == Pi/2, w'[0] == 1}, w, {x, 0, 4 Pi}]
];
Plot[w[x] /. sol, {x, 0, 4 Pi}, AxesOrigin -> {0, 0}]
The equivalent ODE (as given by @Nasser) can be solved by NDSolveValue
when using the {"EquationSimplification"->"Residual"}
Method
:
sol = NDSolveValue[
{
x w''[x] + w'[x] == Sin[w[x]],
w[0] == Pi/2,
w'[0] == 1
},
w,
{x, 0, 12},
Method->{"EquationSimplification"->"Residual"}
];
Plot[sol[x], {x, 0, 12}, AxesOrigin -> {0,0}]