Resolviendo una ODE no lineal de segundo orden con una singularidad en x = 0

Question

Estoy haciendo algunos estudios electromagnéticos nanoestructuras y que tengo que resolver esta ecuación diferencial (exactamente los valores de las constantes no importa, sólo quiero todas las posibles soluciones de y(x) dado que algunos de los valores de estas constantes).

$$ \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) $$

desde x=0 hasta x=R, con las condiciones de frontera

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

Creo que usted no puede encontrar una solución analítica de esta ecuación, por lo que he estado tratando de utilizar métodos numéricos como el método del disparo (dadas las condiciones de contorno, me pareció apropiado).

La cosa es que la singularidad en x=0 no me deja encontrar las soluciones. Puedo obtener diferentes resultados dependiendo de cuántos pasos puedo tomar en el método.

También he publicado esto en la Ciencia Computacional y Física StackExchange, pero por ahora no pude solucionarlo.

Answer 1

El siguiente maxima script le da las fórmulas para el cálculo de los coeficientes de $y_2,\ldots,y_4$ en el ansatz $y=y_1 x + y_2 x^2 + y_3 x^3 + y_4 x^4$. $y_1$ es el parámetro libre para que coincida con el estado inicial en $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]));

El resultado es:

\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*}

Puede utilizar la expansión de la serie de $y(x)$ a huir de la singularidad y dejar que el solucionador numérico patada en algunos suficientemente grande $x$.

También se puede jugar con el orden de la expansión de Taylor para obtener más soluciones exactas o para estimar la calidad de la solución para el escogido $x$.

El cálculo de $\frac{\partial y(x)}{\partial y_1}$ para el compuation de la matriz de un sistema de variacional del sistema que usted necesita para el total-método del disparo debe ser ningún problema.

La expansión de la orden es un parámetro ajustable en la siguiente variante de la anterior secuencia de comandos de maxima. Este se degrada un poco la claridad de la secuencia de comandos, pero se simplifica enormemente jugando con la expansión de la orden.

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

¿Has probado linealización? Si $y$ es pequeña, podemos hacer lo siguiente:

Vamos $\alpha = \frac{K}{A}$, $\beta = \frac{D}{A}$, $\gamma = \frac{\mu HM}{2A}$ y $\delta = -\frac{D}{2A}$ uso de $\sin(2y) = 2\sin(y)\cos(y)$ toda la ecuación puede escribirse como

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

Luego de la reescritura como 2 dimensiones del sistema de orden 1, se obtiene:

$$ \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 $$

Por lo tanto hemos reescrito el sistema no-lineal de la DAE $E(x)y' = A(x)y$$\det(E) = -\frac{1}{x}$. Así que aquí usted puede invertir $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 $$

Y ahora puede pasar esta ODA a sus problemas, y en realidad debería incluso ser analíticamente sovleable.