Raíces escalares y sistemas no lineales

Resolver $f(x)=0$ es una operación central en dinámica. Los ceros de un campo escalar son equilibrios; los ceros de $F\colon\mathbb R^n\to\mathbb R^n$ localizan estados estacionarios; una condición de evento se convierte en una raíz en el tiempo; una órbita periódica puede buscarse como raíz de $P(x)-x$. El algoritmo debe aprovechar la estructura disponible y devolver algo más que un número: residuo, intervalo o paso final, razón de terminación y evidencia de que no se perdió otra solución relevante.

Strogatz y Fuchs formulan la búsqueda de equilibrios y bifurcaciones como problemas de raíces; la implementación se contrasta con la documentación oficial de SciPy. Strogatz, 2015 Steven H. Strogatz (2015) Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering 2.ª ed. · Westview Press Ubicación consultada: caps. 2–3, pp. 15–94

Definición del problema de raíces

Define dominio, unidades y regularidad. Si $f$ tiene una discontinuidad, un cambio de signo puede corresponder a un polo y no a una raíz. Por ejemplo, $f(x)=1/x$ cambia de signo alrededor de cero, pero no está definida allí. Si la evaluación contiene una raíz cuadrada o un logaritmo, el algoritmo debe respetar su dominio.

Un gráfico grueso ayuda a proponer intervalos, pero no certifica unicidad. La derivada puede mostrar monotonicidad: si $f'$ conserva signo y $f(a)f(b)<0$, hay exactamente una raíz en $[a,b]$. Sin ese argumento, un intervalo puede contener varias raíces y la bisección devolverá una de ellas.

El residuo $|f(\widehat x)|$ tampoco equivale al error $|\widehat x-x_\star|$. Cerca de una raíz simple,

$$f(x)\approx f'(x_\star)(x-x_\star),$$

de modo que $|x-x_\star|\approx|f(x)|/|f'(x_\star)|$. Si la derivada es pequeña, un residuo modesto puede ocultar un error grande; si es enorme, un error pequeño puede producir un residuo grande.

Bisección y cota del error

Supón $f$ continua en $[a,b]$ y $f(a)f(b)<0$. El teorema del valor intermedio garantiza al menos una raíz. Se toma $m=(a+b)/2$ y se conserva el subintervalo donde persiste el cambio de signo. Después de $n$ iteraciones, la longitud es

$$\frac{b-a}{2^n}$$

y el punto medio tiene error a lo sumo la mitad de esa longitud respecto de alguna raíz contenida. Para garantizar una semilongitud menor que xtol, basta

$$n\ge \left\lceil\log_2\frac{b-a}{2\,\mathrm{xtol}}\right\rceil.$$

La convergencia es lineal, pero la cota no depende de la derivada. El método falla si no hay cambio de signo, aun cuando exista una raíz de multiplicidad par, y puede encontrar una discontinuidad si no se comprobó el dominio.

Como ejemplo, $f(x)=x^3-x-2$ cambia de signo entre $1$ y $2$. Para una semilongitud menor que $10^{-6}$ se requieren como máximo

$$\left\lceil\log_2(1/(2\cdot10^{-6}))\right\rceil=19$$

iteraciones. El resultado debe incluir el intervalo final; esa pareja de extremos es una evidencia más fuerte que imprimir muchos decimales del punto medio.

Una tabla reproducible conserva en cada iteración $a_n$, $b_n$, $m_n$ y $f(m_n)$. Así se comprueban tanto el cambio de signo como la reducción exacta del intervalo a la mitad.

Iteración de punto fijo y criterio de contracción

La ecuación $f(x)=0$ puede reescribirse como $x=g(x)$ y resolverse mediante

$$x_{k+1}=g(x_k).$$

La equivalencia algebraica no garantiza convergencia. Si existe un intervalo cerrado $I$ tal que $g(I)\subseteq I$ y $|g'(x)|\le q<1$ para todo $x\in I$, el teorema de la contracción asegura un único punto fijo $x_\star$ en $I$. El teorema del valor medio también da la cota

$$|x_{k+1}-x_\star| =|g(x_k)-g(x_\star)| \le q|x_k-x_\star|,$$

por lo que el error decrece al menos geométricamente. Esta justificación muestra por qué dos despejes equivalentes pueden tener comportamientos numéricos opuestos: la derivada de la función iterada, no sólo la ecuación original, determina la estabilidad local del proceso.

El procedimiento reproducible consta de cuatro comprobaciones:

1. elegir $g$ y un intervalo $I$ que contenga la raíz buscada;
2. verificar $g(I)\subseteq I$ y estimar $q=\max_{x\in I}|g'(x)|$;
3. registrar en cada paso $x_k$, $x_{k+1}-x_k$ y el residuo $f(x_{k+1})$;
4. detenerse sólo cuando paso y residuo satisfagan las tolerancias, o declarar que no se verificó la contracción.

Para $f(x)=\cos x-x$ en $I=[0,1]$, la elección $g(x)=\cos x$ cumple $g(I)=[\cos 1,1]\subset I$ y $|g'(x)|=|\sin x|\le\sin 1\approx0.84147$. Con $x_0=0.5$ se obtiene

$$x_1=0.87758256,\quad x_2=0.63901249,\quad x_3=0.80268510,\quad x_4=0.69477803,\quad x_5=0.76819583,$$

y la sucesión converge a $x_\star\approx0.73908513$. Para comparar métodos sobre el mismo problema, bisección conserva el intervalo certificado $[0,1]$, punto fijo conserva la cota de contracción y Newton usa $x_{k+1}=x_k-(\cos x_k-x_k)/(-\sin x_k-1)$. La tabla comparativa debe contar evaluaciones, mostrar el residuo final y distinguir la garantía global del intervalo de la rapidez local.

Velocidad y garantías de Newton y secante

Newton linealiza en $x_k$:

$$f(x_k)+f'(x_k)(x_{k+1}-x_k)=0, \qquad x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}.$$

Cerca de una raíz simple y con suficiente regularidad, el error cumple aproximadamente

$$|e_{k+1}|\approx \left|\frac{f''(x_\star)}{2f'(x_\star)}\right||e_k|^2,$$

lo que da convergencia cuadrática. Lejos de la raíz puede divergir, saltar fuera del dominio o entrar en un ciclo. Una derivada casi nula produce un paso enorme.

La secante reemplaza $f'$ por la pendiente entre dos iterados:

$$x_{k+1}=x_k-f(x_k) \frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}.$$

No requiere derivada y converge superlinealmente cerca de una raíz simple, pero pierde la garantía de un intervalo. Los métodos híbridos mantienen un intervalo con cambio de signo y aceptan pasos rápidos sólo cuando son seguros. La documentación de scipy.optimize.root_scalar expone bisección, Brent, secante, Newton y Halley con contratos distintos. Revisa converged, flag, iterations y function_calls; no uses únicamente .root.

Método de Newton ante raíces múltiples

Para $f(x)=(x-1)^2$, Newton da

$$x_{k+1}=x_k- \frac{(x_k-1)^2}{2(x_k-1)} =\frac{x_k+1}{2}.$$

El error se reduce por un factor $1/2$: la convergencia es lineal, no cuadrática. El cambio de signo desaparece, así que bisección no detecta la raíz en un intervalo genérico. Si se conoce la multiplicidad $m$, la modificación

$$x_{k+1}=x_k-m\frac{f(x_k)}{f'(x_k)}$$

recupera el comportamiento rápido; en la práctica, la multiplicidad suele ser parte del problema que se intenta descubrir.

En bifurcaciones, una raíz doble señala precisamente pérdida de hiperbolicidad. Para $f(x,\mu)=\mu-x^2$, en $\mu=0$ coinciden $f=0$ y $f_x=0$. Fuchs muestra cómo esta degeneración da lugar a la silla-nodo; el deterioro numérico no es accidental, refleja la geometría del problema. Fuchs, 2013 Armin Fuchs (2013) Nonlinear Dynamics in Complex Systems: Theory and Applications for the Life-, Neuro- and Natural Sciences Springer Ubicación consultada: cap. 2, secciones 2.1–2.3, pp. 15–32

Criterios combinados de parada

Un algoritmo puede detenerse cuando el paso satisface

$$|x_{k+1}-x_k|\le \mathrm{xtol}+\mathrm{rtol}|x_{k+1}|,$$

cuando el residuo cae bajo una escala ftol, cuando el intervalo con cambio de signo es pequeño o cuando alcanza un límite de evaluaciones. Conviene combinar paso y residuo. Un paso pequeño puede deberse a estancamiento; un residuo pequeño puede ser insuficiente en un problema mal condicionado.

Escala la función si sus componentes tienen unidades distintas. En un sistema $F(x)=0$, usa una norma de residuo normalizada y una norma de paso. Reporta ambas y el número de iteraciones. Si el método termina por máximo de trabajo, el último iterado no debe etiquetarse como solución.

Método de Newton para sistemas no lineales

Para $F\colon\mathbb R^n\to\mathbb R^n$, Newton resuelve en cada paso

$$J_F(x_k)s_k=-F(x_k), \qquad x_{k+1}=x_k+s_k.$$

No se forma $J^{-1}$: se resuelve el sistema lineal. La condición de $J_F$ influye en la sensibilidad del paso. Cerca de una bifurcación, un jacobiano singular o casi singular puede producir pasos grandes; una estrategia de amortiguamiento usa $x_{k+1}=x_k+\alpha s_k$ con $0<\alpha\le1$ elegido para reducir una función de mérito.

Considere los equilibrios del sistema biestable

$$F(x,v)=\begin{pmatrix}v\\x-x^3-0.45v\end{pmatrix}.$$

El jacobiano es

$$J_F=\begin{pmatrix}0&1\\1-3x^2&-0.45\end{pmatrix}.$$

Semillas cercanas a $(-1,0)$, $(0,0)$ y $(1,0)$ convergen a raíces distintas. Una sola ejecución no enumera todos los equilibrios. Para explorar, usa una malla de semillas, agrupa soluciones por distancia, verifica residuos y después clasifica cada equilibrio con los autovalores del jacobiano dinámico. Hirsch, Smale y Devaney conectan estos ceros con estabilidad y bifurcaciones. Hirsch, 2013 Morris W. Hirsch, Stephen Smale y Robert L. Devaney (2013) Differential Equations, Dynamical Systems, and an Introduction to Chaos 3.ª ed. · Academic Press Ubicación consultada: caps. 1 y 8, pp. 1–20 y 159–186

SciPy ofrece scipy.optimize.root para sistemas, con varios métodos y un objeto de resultado. La rutina no prueba que se hayan encontrado todas las raíces ni que una raíz sea estable; resuelve el problema local planteado desde la semilla.

Registro y validación de la raíz calculada

Una raíz numérica resulta defendible cuando el dominio está declarado, el algoritmo usa la información adecuada, el criterio de parada combina paso y residuo y la búsqueda explica qué otras soluciones podrían haber quedado fuera.