Solving a second-order nonlinear ODE with a singularity on x=0

I'm doing some reasearch on electromagnetic nanostructures and I have to solve this differential equation (the exact values of the constants don't matter, I just want all the possible solutions of y(x) given some values to these constants).

$$\frac{d^2 y}{dx^2}=-\frac{1}{x}\frac{dy}{dx}+\frac{\sin(2y)}{2} (\frac{1}{x^2}+\frac{K}{A})-\frac{D}{A}\frac{\sin(y)}{x}+\frac{\mu HM}{2A}\sin(y)$$

from x=0 till x=R, with the boundary conditions

$$y(0)=0,\ \frac{dy}{dx}(R)=\frac{-D}{2A}$$

I believe you can not find an analitic solution to this equation, so I've been trying to use numerical methods like the shooting method (given the boundary conditions, I found it appropiate).

The thing is that the singularity on x=0 doesn't let me find the solutions. I obtain different results depending on how many steps I take in the method.

I also posted this on Computational Science and Physics StackExchange, but for now I couldn't fix it.

• The condition $y(0)=0$ is (almost) something that follows directly from the differential equation. $x^2 y'' = -xy' + \frac12\sin(2y) + x^2 \frac{K}{2A}\sin(2y) - x\frac{D}{A}\sin(y) + x^2 \frac{\mu HM}{2A}\sin(y)$. For $x=0$ this reduces to $\frac12\sin(2y) = 0$, i.e., $y(0) = \frac{k\pi}2$ with $k\in\mathbb{Z}$. Nov 12, 2016 at 16:06
• Yes, but this changes nothing right? Nov 12, 2016 at 18:51
• It shows the direction. I propose an ansatz for $y$ as Taylor series expansion for small $x$. Use the differential equation in the form from my comment and tune the coefficients of the Taylor series. There should remain one unknown in the ansatz that is required for fitting the other boundary condition via the shooting method. This technique should help you to get sufficiently far away from the singularity. Whether you can work with the Taylor polynomial alone or you need additional numerical integration depends probably on the size of $R$. Nov 12, 2016 at 19:56
• Hey! I've tried it with some solvers, but the solution is quite complicated: wolframalpha.com/input/?i=x%5E2*y%27%27%2Bx*y%27%3Dy(1%2B(a%2Bg)x%5E2-bx) Laguerre and Hypergeometric functions appear, which will take quite long to calculate numerically (I need to do it for all the points close to 0). I can't solve this reduced differential equation through numerical methods either because we didn't get rid of the singularity. At least this method let me continue with the research, but isn't there anything else simpler? Trying to solve this kind of diff. eq. should be common. Nov 13, 2016 at 17:09
• You should write the equation in the form given by me in my previous comment. Otherwise you cannot evaluate the equation at $x=0$. Nov 13, 2016 at 22:00

The following maxima script gives you the formulas for the calculation of the coefficients $y_2,\ldots,y_4$ in the ansatz $y=y_1 x + y_2 x^2 + y_3 x^3 + y_4 x^4$. $y_1$ is the free parameter for matching the initial condition at $x=R$.

declare(x,mainvar);

y : y1*x + y2*x^2 + y3*x^3 + y4*x^4 + y5*x^5;
fDiff : x^2 * diff(y,x,2) + x*diff(y,x,1);
fRHS : sin(2*y)/2 + x^2*K/2/A*sin(2*y) - x*D/A*sin(y) + x^2*mu*H*M/(2*A)*sin(y);

eq : factorout(taylor(fRHS,x,0,4) - fDiff,x);

coeff(eq,x,1);
eq1 : A*scsimp(solve(coeff(eq,x,2));
eq2 : 6*A*scsimp(coeff(eq,x,3));
eq3 : 6*A*scsimp(coeff(eq,x,4));

tex(solve([eq1,eq2,eq3],[y2,y3,y4]));


The result is:

\begin{align*} {\it y_2}&=-{{D\,{\it y_1}}\over{3\,A}}\\ {\it y_3}&=-{{4\,A^2\,{\it y_1}^3+\left(A\,\left(-3\,H\,M\,\mu-6\,K\right)-2\,D^ 2\right)\,{\it y_1}}\over{48\,A^2}}\\ {\it y_4}&={{44\,A^2\,D\, {\it y_1}^3+\left(A\,D\,\left(-11\,H\,M\,\mu-22\,K\right)-2\,D^3 \right)\,{\it y_1}}\over{720\,A^3}} \end{align*}

You can use the series expansion for $y(x)$ to flee from the singularity and let the numerical solver kick in at some sufficiently large $x$.

You can also play with the order of the Taylor expansion to get more exact solutions or to estimate the quality of the solution for the chosen $x$.

Calculating $\frac{\partial y(x)}{\partial y_1}$ for the compuation of the system matrix of the variational system which you need for the overall-shooting method should be no problem.

The expansion order is an adjustable parameter in the following variant of the above maxima script. This degrades a bit the clarity of the script but it greatly simplifies playing around with the expansion order.

declare(x,mainvar);

order : 4;

y : sum(concat('y,i)*x^i,i,1,order);
fDiff : ratsimp(x^2 * diff(y,x,2) + x*diff(y,x,1));
fRHS : sin(2*y)/2 + x^2*K/2/A*sin(2*y) - x*D/A*sin(y) + x^2*mu*H*M/(2*A)*sin(y);

eq : factorout(taylor(fRHS,x,0,order) - fDiff,x);

for i:1 thru order do concat('eq,i) :: coeff(eq,x,i);

sol : solve(makelist(concat('eq,i),i,2,order),makelist(concat('y,i),i,2,order));

tex(sol);


Have you tried linearization? If $y$ is small we can do the following:

Let $\alpha = \frac{K}{A}$, $\beta = \frac{D}{A}$, $\gamma = \frac{\mu HM}{2A}$ and $\delta = -\frac{D}{2A}$ then using $\sin(2y) = 2\sin(y)\cos(y)$ the whole equation can be rewritten as

$$xy' + x^2 y'' = \big((1+\alpha x^2)\underbrace{\cos(y)}_{\approx 1} - \beta x + \gamma x^2 \big)\underbrace{\sin(y)}_{\approx y} \simeq \underbrace{\big( (\alpha + \gamma) x^2 - \beta x +1\big)}_{:=p(x)}y$$

Then rewriting as a 2 dimensional system of order 1 yields:

$$\begin{bmatrix} 0 & 1 \\ x & x^2 \end{bmatrix} \cdot \frac{d}{dx} \begin{bmatrix} y \\ z \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ p(x) & 0 \end{bmatrix} \cdot \begin{bmatrix} y \\ z \end{bmatrix} ,\qquad y(0) = 0, z(R) = \delta$$

Thus we have rewritten the system as a non-linear DAE $E(x)y' = A(x)y$ with $\det(E) = -\frac{1}{x}$. So here you can invert $E$:

$$\frac{d}{dx} \begin{bmatrix} y \\ z \end{bmatrix} = \begin{bmatrix} p(x) & -x \\ 0 & \frac{1}{x} \end{bmatrix} \cdot \begin{bmatrix} y \\ z \end{bmatrix} ,\qquad y(0) = 0, z(R) = \delta$$

And now you can pass this ODE to your solver, and actually it should even be analytically sovleable.

• It looks like something! I'll try to develop a bit more on this way tomorrow, but I'm not sure if y should be small when x -> 0 (maybe if we imply that the second derivative has to be finite at x=0). I'll come back when I have some results, thanks! Nov 12, 2016 at 18:59
• For $x\to 0$ the solution $y(x)$ will definitely be small since you gave the boundary condition $y(0)=0$. That's the reason why the linearization is valid around $x=0$ in the first place. Nov 12, 2016 at 19:23
• True, true, my bad Nov 12, 2016 at 20:20
• Hey! I've tried it with some solvers, but the solution is quite complicated: wolframalpha.com/input/?i=x%5E2*y%27%27%2Bx*y%27%3Dy(1%2B(a%2Bg)x%5E2-bx) Laguerre and Hypergeometric functions appear, which will take quite long to calculate numerically (I need to do it for all the points close to 0). I can't solve this reduced differential equation through numerical methods either because we didn't get rid of the singularity. At least this method let me continue with the research, but isn't there anything else simpler? Trying to solve this kind of diff. eq. should be common. Nov 13, 2016 at 17:09
• The 2nd equation is $z' = \frac{1}{x}z$ which with the boundary condition $z(R) = \delta$ gives simply $z(x) = \frac{\delta}{R}x$. Thus the 1st equation reads $y' = p(x)y - \frac{\delta}{R}$ which is a linear ODE and can easily be solved via the method of an integrating factor. Nov 13, 2016 at 22:12