# Numerically solving an ODE whose right-hand side involves an integral

In a 1955 paper on surfaces of constant negative Gaussian curvature, Marc-Henri Amsler defines a function $w(x)$ by giving the ordinary differential equation $$w'(x)=\frac{1}{x}\int_0^x \sin(w(t))dt,$$ with the initial condition $w(0)=\pi/2$. So the slope $w'(x)$ is the average value of $\sin(w(x))$ over the interval $[0..x]$. In particular, we have $w'(0)=\sin(\pi/2)=1$.

I asked Mathematica to compute $w(x)$ for me numerically by using NDSolve:

NDSolve[{w'[x] == NIntegrate[Sin[w[t]], {t, 0, x}]/x, w[0] == Pi/2}, w[x],
{x, -2, 2}]


but Mathematica (version 10.0.2) responds with a slew of errors, the first batch of which repeatedly report that

NIntegrate::nlim: t = x is not a valid limit of integration. >>


Is this ODE outside of Mathematica's region of expertise, or I am doing something stupid in how I have phrased my Mathematica request?

• That $w(x)$ is a Painlevé transcendent in disguise, which comes from solving the sine-Gordon equation. The difficulties encountered in the two answers below are well-known. Mar 31, 2018 at 19:07

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