Método de Numerov

El método de Numerov (también llamado método de Cowell ) es un método numérico para resolver ecuaciones diferenciales ordinarias de segundo orden en las que no aparece el término de primer orden. Es un método lineal multipaso de cuarto orden . El método es implícito, pero puede hacerse explícito si la ecuación diferencial es lineal.

El método de Numerov fue desarrollado por el astrónomo ruso Boris Vasil'evich Numerov .

El método de Numerov se puede utilizar para resolver ecuaciones diferenciales de la forma

${\displaystyle {\frac {d^{2}y}{dx^{2}}}=-g(x)y(x)+s(x).}$

En él, tres valores de${\displaystyle y_{n-1},y_{n},y_{n+1}}$tomadas en tres puntos equidistantes${\displaystyle x_{n-1},x_{n},x_{n+1}}$se relacionan de la siguiente manera:

${\displaystyle y_{n+1}\left(1+{\frac {h^{2}}{12}}g_{n+1}\right)=2y_{n}\left(1-{\frac {5h^{2}}{12}}g_{n}\right)-y_{n-1}\left(1+{\frac {h^{2}}{12}}g_{n-1}\right)+{\frac {h^{2}}{12}}(s_{n+1}+10s_{n}+s_{n-1})+{\mathcal {O}}(h^{6}),}$

dónde${\ Displaystyle y_ {n} = y (x_ {n})}$,${\displaystyle g_{n}=g(x_{n})}$,${\displaystyle s_{n}=s(x_{n})}$, y${\displaystyle h=x_{n+1}-x_{n}}$.

En física numérica, el método se utiliza para encontrar soluciones de la ecuación de Schrödinger unidimensional para potenciales arbitrarios. Un ejemplo de ello es la resolución de la ecuación radial para un potencial con simetría esférica. En este ejemplo, después de separar las variables y resolver analíticamente la ecuación angular, nos queda la siguiente ecuación de la función radial.${\displaystyle R(r)}$:

${\displaystyle {\frac {d}{dr}}\left(r^{2}{\frac {dR}{dr}}\right)-{\frac {2mr^{2}}{\hbar ^{2}}}(V(r)-E)R(r)=l(l+1)R(r).}$

Esta ecuación se puede reducir a la forma necesaria para la aplicación del método de Numerov con la siguiente sustitución:

${\displaystyle u(r)=rR(r)\Rightarrow R(r)={\frac {u(r)}{r}},}$

${\displaystyle {\frac {dR}{dr}}={\frac {1}{r}}{\frac {du}{dr}}-{\frac {u(r)}{r^{2}}}={\frac {1}{r^{2}}}\left(r{\frac {du}{dr}}-u(r)\right)\Rightarrow {\frac {d}{dr}}\left(r^{2}{\frac {dR}{dr}}\right)={\frac {du}{dr}}+r{\frac {d^{2}u}{dr^{2}}}-{\frac {du}{dr}}=r{\frac {d^{2}u}{dr^{2}}}.}$

Y cuando hacemos la sustitución, la ecuación radial se convierte en:

${\displaystyle r{\frac {d^{2}u}{dr^{2}}}-{\frac {2mr}{\hbar ^{2}}}(V(r)-E)u(r)={\frac {l(l+1)}{r}}u(r),}$

o

${\displaystyle -{\frac {\hbar ^{2}}{2m}}{\frac {d^{2}u}{dr^{2}}}+\left(V(r)+{\frac {\hbar ^{2}}{2m}}{\frac {l(l+1)}{r^{2}}}\right)u(r)=Eu(r),}$

lo cual es equivalente a la ecuación de Schrödinger unidimensional, pero con el potencial efectivo modificado.

${\displaystyle V_{\text{eff}}(r)=V(r)+{\frac {\hbar ^{2}}{2m}}{\frac {l(l+1)}{r^{2}}}=V(r)+{\frac {L^{2}}{2mr^{2}}},\quad L^{2}=l(l+1)\hbar ^{2}.}$

Podemos proceder a resolver esta ecuación del mismo modo que habríamos resuelto la ecuación de Schrödinger unidimensional. Podemos reescribir la ecuación de forma ligeramente diferente y así ver con mayor claridad la posible aplicación del método de Numerov:

${\displaystyle {\frac {d^{2}u}{dr^{2}}}=-{\frac {2m}{\hbar ^{2}}}(E-V_{\text{eff}}(r))u(r),}$

${\displaystyle g(r)={\frac {2m}{\hbar ^{2}}}(E-V_{\text{eff}}(r)),}$

${\displaystyle s(r)=0.}$

Se nos da la ecuación diferencial

${\displaystyle y''(x)=-g(x)y(x)+s(x).}$

Para derivar el método de Numerov para resolver esta ecuación, comenzamos con el desarrollo de Taylor de la función que queremos resolver,${\displaystyle y(x)}$, alrededor del punto${\displaystyle x_{0}}$:

${\displaystyle y(x)=y(x_{0})+(x-x_{0})y'(x_{0})+{\frac {(x-x_{0})^{2}}{2!}}y''(x_{0})+{\frac {(x-x_{0})^{3}}{3!}}y'''(x_{0})+{\frac {(x-x_{0})^{4}}{4!}}y''''(x_{0})+{\frac {(x-x_{0})^{5}}{5!}}y'''''(x_{0})+{\mathcal {O}}(h^{6}).}$

Denotando la distancia desde${\displaystyle x}$a${\displaystyle x_{0}}$por${\displaystyle h=x-x_{0}}$, podemos escribir la ecuación anterior como

${\displaystyle y(x_{0}+h)=y(x_{0})+hy'(x_{0})+{\frac {h^{2}}{2!}}y''(x_{0})+{\frac {h^{3}}{3!}}y'''(x_{0})+{\frac {h^{4}}{4!}}y''''(x_{0})+{\frac {h^{5}}{5!}}y'''''(x_{0})+{\mathcal {O}}(h^{6}).}$

Si discretizamos el espacio de manera uniforme, obtenemos una cuadrícula de${\displaystyle x}$puntos, donde${\displaystyle h=x_{n+1}-x_{n}}$. Al aplicar las ecuaciones anteriores a este espacio discreto, obtenemos una relación entre el${\displaystyle y_{n}}$y${\displaystyle y_{n+1}}$:

${\displaystyle y_{n+1}=y_{n}+hy'(x_{n})+{\frac {h^{2}}{2!}}y''(x_{n})+{\frac {h^{3}}{3!}}y'''(x_{n})+{\frac {h^{4}}{4!}}y''''(x_{n})+{\frac {h^{5}}{5!}}y'''''(x_{n})+{\mathcal {O}}(h^{6}).}$

Desde el punto de vista computacional, esto equivale a dar un paso adelante en una cantidad${\displaystyle h}$. Si queremos dar un paso atrás , reemplazamos cada${\displaystyle h}$con${\displaystyle -h}$y obtener la expresión para${\displaystyle y_{n-1}}$:

${\displaystyle y_{n-1}=y_{n}-hy'(x_{n})+{\frac {h^{2}}{2!}}y''(x_{n})-{\frac {h^{3}}{3!}}y'''(x_{n})+{\frac {h^{4}}{4!}}y''''(x_{n})-{\frac {h^{5}}{5!}}y'''''(x_{n})+{\mathcal {O}}(h^{6}).}$

Tenga en cuenta que solo las potencias impares de${\displaystyle h}$experimentó un cambio de signo. Sumando las dos ecuaciones, obtenemos que

${\displaystyle y_{n+1}-2y_{n}+y_{n-1}=h^{2}y''_{n}+{\frac {h^{4}}{12}}y''''_{n}+{\mathcal {O}}(h^{6}).}$

Podemos resolver esta ecuación para${\displaystyle y_{n+1}}$sustituyendo la expresión dada al principio, es decir${\displaystyle y''_{n}=-g_{n}y_{n}+s_{n}}$. Para obtener una expresión para el${\displaystyle y''''_{n}}$factor, simplemente tenemos que diferenciar${\displaystyle y''_{n}=-g_{n}y_{n}+s_{n}}$dos veces y aproximámoslo de nuevo de la misma manera que lo hicimos anteriormente:

${\displaystyle y''''_{n}={\frac {d^{2}}{dx^{2}}}(-g_{n}y_{n}+s_{n}),}$

${\displaystyle h^{2}y''''_{n}=-g_{n+1}y_{n+1}+s_{n+1}+2g_{n}y_{n}-2s_{n}-g_{n-1}y_{n-1}+s_{n-1}+{\mathcal {O}}(h^{4}).}$

Si ahora sustituimos esto en la ecuación anterior, obtenemos

${\displaystyle y_{n+1}-2y_{n}+y_{n-1}={h^{2}}(-g_{n}y_{n}+s_{n})+{\frac {h^{2}}{12}}(-g_{n+1}y_{n+1}+s_{n+1}+2g_{n}y_{n}-2s_{n}-g_{n-1}y_{n-1}+s_{n-1})+{\mathcal {O}}(h^{6}),}$

o

${\displaystyle y_{n+1}\left(1+{\frac {h^{2}}{12}}g_{n+1}\right)-2y_{n}\left(1-{\frac {5h^{2}}{12}}g_{n}\right)+y_{n-1}\left(1+{\frac {h^{2}}{12}}g_{n-1}\right)={\frac {h^{2}}{12}}(s_{n+1}+10s_{n}+s_{n-1})+{\mathcal {O}}(h^{6}).}$

Esto da como resultado el método de Numerov si ignoramos el término de orden.${\displaystyle h^{6}}$De ello se deduce que el orden de convergencia (suponiendo estabilidad) es 4.

• Hairer, Ernst; Norsett, Syvert Paul; Wanner, Gerhard (1993), Resolución de ecuaciones diferenciales ordinarias I: problemas no rígidos , Berlín, Nueva York: Springer-Verlag , ISBN 978-3-540-56670-0Este libro incluye las siguientes referencias :
• Numerov, Boris Vasil'evich (1924), "Un método de extrapolación de perturbaciones", Monthly Notices of the Royal Astronomical Society , 84 (8): 592– 601, Bibcode : 1924MNRAS..84..592N , doi : 10.1093/mnras/84.8.592.
• Numerov, Boris Vasil'evich (1927), "Nota sobre la integración numérica de d ² x /d t ² = f ( x , t )", Astronomische Nachrichten , 230 (19): 359– 364, Bibcode : 1927AN....230..359N , doi : 10.1002/asna.19272301903.