8
$\begingroup$

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?

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Mar 31, 2018 at 19:07

3 Answers 3

8
$\begingroup$

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

Mathematica graphics

It differs from Nasser's answer by less than 3 * 10^-8 over the interval.

$\endgroup$
6
$\begingroup$

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

Mathematica graphics

$\endgroup$
4
$\begingroup$

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.