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

Question

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.

Answer 1

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);

Answer 2

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.