Todas las entradas de: cj817290-ovh

Problemas de probabilidades

1. En una red local hay las conexiones mostradas en la figura, donde los números indican las probabilidades de que cada rama esté abierta en un cierto intervalo de tiempo dado. Suponiendo que las probabilidades son independientes entre sí, calcular la probabilidad de que haya transmisión de datos entre A y D por cualquier camino. Suponiendo que hay transmisión de datos entre A y D, calcular la probabilidad de que se esté transmitiendo por la ruta ACD.

graf_probabilitats

 Para transmitir entre A y D hay tres caminos: ABD, AD, ACD; si cualquiera de ellos está abierto, hay comunicación entre A y D. ¿cuáles son las posibilidades?

  • ABD abierto, AD y ACD cerrados
  • AD abierto, ABD y ACD cerrados
  • ACD abierto, ABD y AD cerrados
  • ABD y AD abiertos, ACB cerrado
  • ...
  • ABD, AD y ACB abiertos

Vemos que hay bastantes posibilidades a considerar; en estos casos es conveniente pensar en el suceso contrario: ¿cuándo no habrá transmisión entre A y D? Sólo cuando ABD, AD y ACB estén todos cerrados. El camino ABD estará cerrado si AB lo está, o bien BD lo está; teniendo en cuenta que los sucesos son independientes, la probabilidad de "ABD cerrado" es:

\begin{array}{l}\text{P}\left(\text{ABD cerrado}\right)\;=\text{P}\left(\text{AC cerrado}\cup\text{CD cerrado}\right)\;=\text{P}\left(\text{AC cerrado}\right)+\text{P}\left(\text{CD cerrado}\right)\text{-P}\left(\text{AC cerrado}\cap\text{CD cerrado}\right)=\\0.1+0.2-0.1\cdot0.2=0.28\end{array}

ya que P(A cerrado) = 1 - P(A abierto) = 1 - 0.9, e idénticamente para B. La probabilidad de "ACD cerrado" es numéricamente la misma:

P(ACD \; cerrado) = P(AC \; cerrado \cup CD \; cerrado) = P(AC \; cerrado) + P(CD \; cerrado) - P(AC \; cerrado \cap  CD \; cerrado) = 0.2 + 0.1 - 0.2 = 0.28.

Entonces P(no se transmite entre A y D)=P(ABD, AD y ACB todos cerrados)  = P(ABD cerrado \cap AD cerrado \cap ACB cerrado) = P(ABD cerrado})·P(AD cerrado)·P(ACB cerrado) = 0.28·0.3·0.28 = 0.02352.

Por tanto P(se transmite entre A y D) = 1 - P(no se transmite entre A y D) = 1 - 0.02352 = 0.97648.

NOTA: puede ser didáctico realizar simulaciones de probabilidades con hoja de cálculo para verificar experimentalmente los cálculos. En este ejercicio es simple de hacer: usando la función aleatorio() que llevan todas las hojas de cálculo, y con la función lógica =SI(condición; valor_si_cierto;valor_si falso), se puede crear una hoja que presente el valor 1 siempre que el valor aleatorio esté en el intervalo [0,p] siendo p las probabilidades dadas de la red:

AB BD AD AC CD
1 1 1 0 1
1 1 1 1 1
0 0 0 0 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
0 1 0 1 1

Así, por ejemplo, la columna AB presenta un 1 siempre que en esa casilla se haya generado un valor aleatorio en el intervalo [0, 0.9]; cuando hay un 1 significa que la ruta AB está abierta, con un 0 está cerrada. Observad que los valores de esta tabla son binarios, con 1=cierto (hay transmisión), 0=falso (no hay transmisión).

Ampliamos ahora con más columnas:

AB y BD AC y CD (AB·BD)+AD+(AC·CD)
1 0 2
1 1 3
0 0 0
1 1 3
1 1 3
1 1 3
1 1 3
1 1 3
0 1 1

En "AB y AD" multiplicamos las columnas AB por AD, en "AC y AD" lo mismo, pero en la columna (AB·BD)+AD+(AC·CD) al sumar no obtenemos un número binario: si éste valor es cero significa que no hay transmisión entre A y D (todo cerrado) y con un valor superior a cero hay transmisión entre A y D (alguna ruta abierta). Contando el número de celdas superiores a cero de ésta última columna y dividiendo por el número de filas obtenemos una estimación de la probabilidad pedida, tanto mejor como más filas haya. Con 1000 filas se obtienen valores del orden de 0.98.

separador2

 2.  Se elige al azar un número de 3 bits x_1x_2x_3 donde x_i=0x_i=1. Definimos las variables aleatorias X: el número de ceros que tienen conjuntamente los dos primeros bits, Y: número total de unos entre los tres bits. Calcular la tabla de distribución conjunta de probabilidad de X,Y. ¿Cuál es la covarianza de X,Y?

 Para calcular las probabilidades conjuntas tenemos que saber primero las posibles combinaciones de valores, son: X puede valer {0, 1, 2}, Y puede valer {0, 1, 2, 3}. Por tanto tendremos 3·4 = 12 combinaciones de valores (X,Y), que son {(0,0), (0,1), (0,2), (0,3), ..., (2,2), (2,3)}.

Vamos por las probabilidades conjuntas: dado un suceso A de la variable X y un suceso B de la variable Y, la probabilidad conjunta de A y B es P\left(A\cap B\right)=P(A\vert B)\cdot P(B), que será igual a P(A)·P(B) sólo si A, B son sucesos independientes; en este caso, no podemos presuponer independencia entre X, Y, luego aplicamos la primera igualdad.

Veamos un ejemplo de cálculo: sea A el suceso X=1, B el suceso Y=2; la probabilidad condicionada P(X = 1 | Y = 2) se obtiene considerando todos los casos Y=2 y viendo en que proporción de ellos se cumple X=1.

Si Y=2 los bits x_1x_2x_3={(1,1,0), (1,0,1), (0,1,1)}, observando los dos primeros bits, vemos que en dos casos, {(1,0,1), (0,1,1)}, tenemos un cero en uno de los bits; luego la proporción de casos X=1 dentro de Y=2 es de 2/3.

Calculamos ahora la P(Y = 2) como la proporción de casos en que tenemos
2 unos respecto al total de combinaciones de los tres bits, que son 2³=8; la proporción es pues 3/8. Por tanto, P\left(A\cap B\right)=\frac23\cdot\frac38=\frac28=\frac14. Por otro lado, la P(X = 1) se obtiene con la proporción de casos en que tenemos 1 cero en los dos primeros bits respecto al total de casos en esos dos primeros bits: {(1,0),(0,1),(1,1),(0,0)}, por tanto es P(X = 1) = 3/4. La probabilidad P(X = 1)·P(Y = 2) = 3/4 · 3/8 = 9/32 que es distinta de la obtenida para P\left(A\cap B\right), luego los sucesos no son independientes.

Obviamente no vamos a hacer un cálculo tan largo para las otras 11 combinaciones (X,Y), lo abreviamos haciendo una discusión de casos simples:

  • Si el número de unos es de 3 (Y=3), entonces no hay ningún cero, luego X debe valer cero con seguridad (X=1 con probabilidad 100%).
  • Si el número de unos es de 0 (Y=0), entonces todos los bits son cero, luego X debe valer 2 con seguridad (X=2 con probabilidad 100%).
  • El caso X=0, Y=0 es imposible (sucesos incompatibles), ya que X=0 implica que hay dos unos en los dos primeros bits, luego Y tiene que valer al menos 2; por tanto el caso X=0, Y=1 también es imposible.
  • Si Y=2 hay un sólo cero en los tres bits, luego X no puede valer 2; tenemos que P(X = 2 | Y = 2) = 0.
  • Si Y=1, los bits han de ser {(1,0,0), (0,1,0), (0,0,1)}; luego P(X = 1|Y = 1) = 2/3, P(X = 2|Y = 1) = 1/3.

Resumimos todo lo que tenemos en una tabla de probabilidades condicionadas P(X | Y), marcamos las casillas que hemos visto que tienen probabilidad 0 (sucesos incompatibles):

Probabilidades condicionadas P(X | Y)
Probabilidades condicionadas P(X | Y)

Para obtener la tabla de probabilidades conjunta P\left(A\cap B\right) usaremos la fórmula P\left(A\cap B\right)=P(A\vert B)\cdot P(B) y algunas propiedades útiles:

  1. la suma de probabilidades por filas (distribución marginal de X) coincide con las probabilidades P(X),
  2. la suma de probabilidades por columnas (distribución marginal de Y) coincide con las probabilidades P(Y), y
  3. la suma total ha de ser 1.

Para aplicarlo, será útil tener la tabla de probabilidades para la variable Y:

Y 0 1 2 3
P(Y) 1/8 3/8 3/8 1/8

y también para la variable X:

X 0 1 2
P(X) ¼ ½ ¼

Obtenemos la tabla conjunta X,Y:

Probabilidades conjuntas
Probabilidades conjuntas

Las casillas en azul se han obtenido aplicando las propiedades 1 y 2, no ha sido necesario el cálculo de probabilidades.

Para calcular la covarianza usamos la fórmula:

Cov\left(X,Y\right)=\sum_x\sum_y\left(x-\mu_x\right)\left(y-\mu_y\right)P\left(x,y\right)

Necesitamos los valores medios de las variables:

\mu_x=\underset x{\sum x}\cdot P(x),\;\mu_y=\underset y{\sum y}\cdot P(y)

Los obtenemos de las tablas de probabilidades para X e Y:

Y 0 1 2 3
P(Y) 1/8 3/8 3/8 1/8
Y·P(Y) 0 0.375 0.75 0.375 suma=1.5
X 0 1 2
P(X) 1/4 1/2 1/4  
X·P(X) 0 0.5 0.5 suma=1

Los elementos que entran en el cálculo de la covarianza los disponemos también en forma de tabla:

Y
X 0 1 2 3
0 0 0 -0.0625 -0.1875
1 0 0 0 0
2 -0.1875 -0.0625 0 0

La suma de todos ellos es la covarianza: -0.5, un valor negativo indica dependencia inversa:  valores de X  grandes implican pequeños valores de Y.

separador2

Diferenciación e integración numérica

Introducción

En las aplicaciones prácticas, a menudo sucede que tenemos valores de la variación de una función en algunos puntos sin tener su expresión analítica, y queremos inferir de ellos la función derivada en algún otro punto. De la misma forma, podemos necesitar la integral definida de esa función de la que sólo sabemos sus valores en algunos puntos; en el caso de la integral además, puede darse el caso de que, incluso teniendo la expresión analítica de la función a integrar, y siendo la función integrable, no exista la función primitiva. En estos casos la solución pasa por aproximar la derivada y la integral usando los métodos del cálculo numérico.

Ejemplo 1: Por experiencia previa sabemos que la velocidad media de transmisión de datos V(t) de una línea de comunicaciones digitales es, en función de la hora, la siguiente:

Hora:

10:00 12:00 14:00 16:00 18:00 20:00

Velocidad (Mbits/s):

80 70 68 69 72 80

Nos interesa conocer (a) el valor de la función derivada V'(t) para el tiempo t=15:00, (b) el valor de la integral definida de la función  V(t) en el intervalo de tiempo [10:00, 20:00]. La derivada V'(t) puede usarse para estudiar el porqué de las variaciones de la velocidad de transmisión durante el día, y la integral está relacionada con el tráfico total de datos a lo largo del día.

Derivación aproximada por interpolación

 Si tenemos una función continua y derivable en un intervalo [a,b] de la cual tenemos sus valores en n+1 puntos diferentes x_0, x_1, \cdots, x_n, y queremos saber el valor de su derivada en algún punto x del intervalo, una forma posible de hacerlo es:

  1. calcular el polinomio interpolador que aproxima a la función en los puntos dados
  2. derivar el polinomio y evaluar la derivada en x

Los detalles del cálculo del polinomio interpolador pueden verse en el post Interpolación polinómica.

Ejemplo 2: evaluar un valor aproximado para la función derivada V'(t) en el tiempo t=15:00 de la función del ejemplo 1. Formamos la tabla de interpolación:

x f(x) Dif 1 Dif 2 Dif 3 Dif 4 Dif 5
10 80
12 70 -5
14 68 -1 -1
16 69 0.5 -0.375 0.1041666667
18 72 1.5 -0.25 0.0208333333 -0.0104166667
20 80 4 -0.625 -0.0625 -0.0104166667 0

El polinomio interpolador es:

\begin{array}{l}80-5(x-10)-1(x-10)(x-12)+0.1041666667(x-10)(x-12)(x-14)+(x-10)(x-12)(x-14)(x-16)\times(-0.0104166667)\\=-0.0104167x^4+0.645833x^3-15.2083x^2+150.417x-445\end{array} y la derivada es -0.0416667x^3+1.9375x^2-30.4167x+150.417,

sustituyendo el valor x=15 obtenemos f'(15)=-10.5.

Polinomio interpolador
Polinomio interpolador

Derivada del polinomio
Derivada del polinomio

Ejemplo 3: De una función de la que desconocemos su expresión analítica sabemos su valor en tres puntos, son: \left(0,1\right),\;\left(\frac{2\mathrm\pi}6,0.5\right),\;\left(\frac{\mathrm\pi}2,0\right). Evaluar el valor de la derivada de la función en el punto x=\frac{\mathrm\pi}4.

Formamos la tabla de interpolación:

 

x f(x) f\left[x_{i-1},x_i\right] f\left[x_{i-2},x_{i-1},x_i\right]
0 1
\frac{2\mathrm\pi}6 0.5 -0.477
\frac{\mathrm\pi}2 0 -0.955 -0.304

El polinomio interpolador es:

P(x)=1-0.477\left(x-0\right)-0.304\left(x-0\right)\left(x-\frac{2\mathrm\pi}6\right)=-0.304x^2-0.159x+1,

derivando y sustituyendo el valor pedido: P'(x)=-0.608x-0.159;\;P'(\frac{\mathrm\pi}4)=-0.608\frac{\mathrm\pi}4-0.159=-0.637;, por tanto el valor aproximado es f'(\frac{\mathrm\pi}4)\approx-0.637;. Los valores del enunciado realmente corresponden a la función y=cos(x), con derivada y'=-sin(x), y el valor exacto es -sin(\frac{\mathrm\pi}4)=-0.707, la precisión obtenida no es muy alta debido a que sólo hemos usado tres puntos de interpolación.

Integración aproximada por interpolación

Dada una función integrable Riemann f(x) de la cual no podemos obtener su función primitiva F\left(x\right)=\int f\left(x\right), pero sí sabemos sus valores en un conjunto de puntos interior al intervalo (a,b), entonces podemos calcular la integral definida F\left(x\right)=\int_a^bf\left(x\right) con el siguiente método aproximado: consideramos una partición del intervalo (a,b) según n+1 puntos equidistantes x_0=a, x_1, ..., x_{n-1}, x_n=b;  podemos aproximar la función f(x) en cada subintervalo (x_{i-1},x_i) por su polinomio interpolador de grado 1, o sea una recta. El área comprendida entre la gráfica de la función y el eje X podrá aproximarse sumando las áreas de los trapecios formados por los puntos x_{i-1},x_i, f(x_{i-1}),f(x_i), como muestra la figura:

integración por interpolación de grado 1: trapecios
integración por interpolación
de grado 1: trapecios. Se muestra la gráfica de f(x) y los segmentos rectos que aproximan a la función en cada intervalo, así como uno de los trapecios: entre los puntos segundo y tercero.

El área del trapecio formado por esos puntos es:

A_{i-1,i}=f\left(x_{i-1}\right)\cdot\triangle x+\frac12\left|f\left(x_{i-1}\right)-f\left(x_i\right)\right|\cdot\triangle x=\frac12\left|f\left(x_{i-1}\right)+f\left(x_i\right)\right|\cdot\triangle x=\frac12\left(f\left(x_{i-1}\right)+f\left(x_i\right)\right)\cdot\triangle x

donde hemos denominado \triangle x a la longitud de los subintervalos. Sumando todos los n trapecios obtenemos la denominada fórmula del trapecio para la integral:

\begin{array}{l}\int_a^bf\left(x\right)\operatorname dx\approx\sum\nolimits_1^nA_{i-1,i}=\sum\nolimits_1^n\frac12\left(f\left(x_{i-1}\right)+f\left(x_i\right)\right)\cdot\triangle x=\\\frac{\triangle x}2\left(f\left(x_0\right)+2f\left(x_1\right)+2f\left(x_2\right)+\dots+2f\left(x_{n-1}\right)+f\left(x_n\right)\right)\end{array}

Ejemplo 4: Para la función f(x)=x^3-x^2+1 aproximamos su integral definida I en el intervalo [0, 1] con n = 4 subintervalos, obtenemos:

valores de f(x) en los puntos de interpolación:

x 0 0,25 0,5 0,75 1
y 1 0,953125 0,875 0,859375 1

anchura de cada subintervalo: (1 - 0)/4 = 0.25; aplicamos la fórmula del trapecio:

\int_0^1f\left(x\right)\operatorname dx\approx\frac{0.25}2\left(1+2\cdot0.953+2\cdot0.875+2\cdot0.859+1\right)=0.922

Si calculamos el valor exacto, usando la función primitiva, hallamos el valor 0.91\overset\frown6.

Fórmula de Simpson

Con el mismo procedimiento de dividir el intervalo de integración en n subintervalos, pero usando polinomios interpoladores de grado 2 en cada subintervalo, obtenemos la denominada fórmula de Simpson, omitimos los detalles del cálculo:

\int_a^bf\left(x\right)\operatorname dx\approx\frac{\triangle X}3\left(f\left(x_0\right)+4f\left(\frac{x_0+x_1}2\right)+2f\left(x_1\right)+4f\left(\frac{x_1+x_2}2\right)+2f\left(x_2\right)+\dots+2f\left(x_{n-1}\right)+4f\left(\frac{x_{n-1}+x_n}2\right)+f\left(x_n\right)\right).

En cada subintervalo se usan tres puntos: los extremos y el punto medio, para n subintervalos tenemos un total de 2n+1 puntos. Lógicamente esperamos que la fórmula de Simpson sea más precisa que la de los trapecios.

Ejemplo 5: Obtener \int_0^1\ln\left(\cos^2\left(x\right)\right)\operatorname dx por el método de Simpson, con n=4.

Calculamos los valores de la función que aparecen en la fórmula de Simpson:

x 0 0,125 0,25 0,375 0,5 0,625 0,75 0,875 1
y 0 -0,0156658605 -0,0631621025 -0,1440500235 -0,2611684809 -0,4190654025 -0,6247997958 -0,8894614471 -1,2312529408

Sustituimos los valores en la fórmula:

\int_0^1\ln\left(\cos^2\left(x\right)\right)\operatorname dx\approx\frac{0.25}6\left[0+4\cdot\left(-0,01567\right)+2\cdot\left(-0,0632\right)+\dots+\left(-1,2313\right)\right]=-0,3751.

Problemas de EDOs de primer orden

Enunciados

1. Resolver la ecuación lineal de primer orden y'x+y=x\sin(x).

2. Resolver la ecuación de Bernouilli y'+\frac yx=y^2\frac{x\cos\left(x\right)-\sin\left(x\right)}x

3.  Integrar la ecuación de Lagrange y=-x+\frac{1-y'}{1+y'}.

4. Integrar la ecuación diferencial exacta \left(2x+\frac1y\right)\operatorname dx+\left(\frac1y-\frac x{y^2}\right)\operatorname dy=0.

 

Soluciones

separador2

 

 

1. Resolver la ecuación lineal de primer orden y'x+y=x\sin(x).

Solución: Ponemos la ecuación en la forma estándard y'+P(x)y=Q(x), dividiendo todo por x: y'+(1/x)y=\sin(x), identificamos P(x)=1/x, Q(x)=\sin(x). Calculamos u=exp(-\int e^{P(x)}\operatorname dx), \;v=\int Q(x)\cdot\left[\int e^{P(x)}\operatorname dx\right]\operatorname dx y la solución será y=u(C+v):

\begin{array}{l}u=e^{\int-\frac1x\operatorname dx}=e^{-\ln\left(x\right)}=\frac1x;\;\\v=\int\sin\left(x\right)\left[e^{\int\frac1x\operatorname dx}\right]\operatorname dx=\int\sin\left(x\right)\cdot x\operatorname dx=\sin\left(x\right)-x\cos\left(x\right)\end{array}

la última integral se puede hacer por partes. La integral general queda:

y=u\left(C+v\right)=\frac1x\left(C+\sin\left(x\right)-x\cos\left(x\right)\right).

separador2

 

 

2. Resolver la ecuación de Bernouilli y'+\frac yx=y^2\frac{x\cos\left(x\right)-\sin\left(x\right)}x

Solución: Con el cambio y=1/u se transforma en una ecuación lineal:

v'-\frac1xv=-\frac{x\cos\left(x\right)-\sin\left(x\right)}x

Para resolver la ecuación lineal, que es del tipo y'+P(x)y=Q(x), hallamos v=-\int e^{P(x)}\operatorname dx,\;w=\int Q(x)\cdot\left[\int e^{P(x)}\operatorname dx\right]\operatorname dx y la solución será y=v(C+w); obtendremos: u=Cx-\sin(x), deshaciendo el cambio, y=(Cx-\sin(x))^{-1}.

separador2

 

 

3.  Integrar la ecuación de Lagrange y=-x+\frac{1-y'}{1+y'}.

Solución: Las ecuaciones de Lagrange por lo general no pueden resolverse de forma explícita, esto es, dando una expresión y=f(x), si no que ha darse en forma paramétrica: x=f(p), y=g(p). El procedimiento es como sigue:

1. Hacemos el cambio y'=p, resulta: y=-x+\frac{1-p}{1+p}=F(x)+G(p)

2. Derivamos toda la ecuación respecto de x, teniendo en cuenta que \frac{\operatorname d{G(p)}}{\operatorname dx}=\frac{\operatorname d{G(p)}}{\operatorname dp}\frac{\operatorname dp}{\operatorname dx}, resulta: p=-1+\frac{-1\left(1+p\right)-\left(1-p\right)1}{\left(1+p\right)^2}\frac{dp}{dx}

 

3. Reordenamos, considerando x como función de p, obtenemos una ecuación lineal en x(p), que en este caso es además separable, que podemos resolver:

y=-x+\frac{1-p}{1+p}=F(x)+G(p)

4. Sustituyendo esta expresión de x en la ecuación original, cambiando y' por p, obtenemos la y en función del parámetro p:

\begin{array}{l}y=-x+\frac{1-p}{1+p}=\frac{-1}{\left(1+p\right)^2}-C+\frac{1-p}{1+p}=-C+\frac{-1+\left(1-p\right)\left(1+p\right)}{\left(1+p\right)^2}=\\-C-\frac{p^2}{\left(1+p\right)^2}\\\\\end{array}

El conjunto de la dos igualdades para x, y nos da el haz de curvas solución en función del parámetro p=y', que es la pendiente de la curva:

\left.\begin{array}{r}x=\frac1{\left(1+p\right)^2}+C\\y=-C-\frac{p^2}{\left(1+p\right)^2}\end{array}\right\}

separador2

 

 

4. Integrar la ecuación diferencial exacta \left(2x+\frac1y\right)\operatorname dx+\left(\frac1y-\frac x{y^2}\right)\operatorname dy=0.

Solución: Comprobamos que es exacta: la ecuación M\left(x,y\right)\operatorname dx+N\left(x,y\right)\operatorname dy=0 es exacta cuando \frac{\partial M}{\partial y}=\frac{\partial N}{\partial x}, en efecto:

\begin{array}{l}\frac{\partial M}{\partial y}=\frac{\partial\left(2x+\frac1y\right)}{\partial y}=-\frac1{y^2};\\\frac{\partial N}{\partial x}=\frac{\partial\left(\frac1y-\frac x{y^2}\right)}{\partial x}=-\frac1{y^2};\end{array}

Integramos los dos miembros de la ecuación, respecto a x o a y:

\begin{array}{l}\int\left(2x+\frac1y\right)\operatorname dx=x^2+\frac xy+C\left(y\right);\\\int\left(\frac1y-\frac x{y^2}\right)\operatorname dy=\ln\left(y\right)+\frac xy+C\left(x\right);\end{array}

Igualamos ambos resultados:

\begin{array}{l}x^2+\frac xy+C\left(y\right)=\ln\left(y\right)+\frac xy+C\left(x\right)\Leftrightarrow\left\{\begin{array}{l}C\left(x\right)=x^2\\C\left(y\right)=\ln\left(y\right)\end{array}\right.\\\end{array}

La solución es F(x,y)=C donde F\left(x,y\right)=\int M\operatorname dx=\int N\operatorname dy=x^2+\frac xy+\ln\left(y\right), o sea, queda en forma implícita, no damos una expresión y=f(x), sino:

x^2+\frac xy+\ln\left(y\right)=C

separador2

 

 

 

 

EDOs lineales de orden superior

separador2

EDOs lineales de segundo orden

Tienen la forma A(x)y''+B(x)y'+C(x)y=F(x). En el caso de que F(x)=0 se llaman ecuaciones homogéneas, que discutimos ahora.

Ecuaciones lineales homogéneas de segundo orden

La ecuación homogénea A(x)y''+B(x)y'+C(x)y=0 asociada a la ecuación lineal completaA(x)y''+B(x)y'+C(x)y=F(x) tiene la siguiente propiedad:

Proposición 1 (principio de superposición): Si y_{1},y_{2} son dos soluciones independientes de la ecuación lineal, entonces y=c_{1}y_{1}+c_{2}y_{2} es la solución general de la ecuación lineal, siendo c_{1},c_{2} valores reales cualesquiera.

Observemos que el principio de superposición se aplica por igual a ecuaciones lineales homogéneas y completas. Además, exige que las dos soluciones sean independientes. Después del siguiente ejemplo lo aclaramos.

Ejemplo 1:  y_{1}=\cos x e y_{2}=\sin x son dos soluciones independientes de la ecuación \mbox{y''+y=0}, por tanto y=c_{1}\cos x+c_{2}\sin x es la solución general.

separador2

Combinaciones lineales de soluciones e independencia lineal

Las funciones y_{1}=\cos x e y_{2}=2\cos x som ámbas soluciones de la ecuación del ejemplo 1, pero no son linealmente independientes, pues podemos escribir la segunda en función de la primera: y=C\sin\left(x\right). En general, si tenemos dos funciones y_1,y_2 que pueden escribirse una en función de la otra, podremos ponerlas en la forma C_1y_1+C_2y_2=0; derivando respecto de x obtenemos una segunda condición: C_1y'_1+C_2y'_2=0. Las dos ecuaciones anteriores forman un sistema de ecuaciones con dos incógnitas:

\left.\begin{array}{r}C_1y_1+C_2y_2=0\\C_1y'_1+C_2y'_2=0\end{array}\right\}

Si este sistema tiene soluciones C_1, C_2, entonces las funciones y_1,y_2 no son independientes; esto sucederá cuando el determinante del sistema no valga cero; por tanto, la condición de independencia será que ese determinante valga cero:

\begin{bmatrix}y_1&y_2\\y'_1&y'_2\end{bmatrix}=0\Rightarrow y_1,y_2\;\text{linealmente independientes}

Al determinante anterior formado con las funciones y_1,y_2 y sus derivadas se le conoce por el nombre de determinante Wronskiano del conjunto de funciones {y_1,y_2}, y se generaliza fácilmente a conjuntos de n funciones \left\{y_1,y_{2,\dots,}y_n\right\}.

Ejemplo 2: Las funciones e^x, e^{-x},Cosh(x) no son linealmente independientes, pues su Wronskiano es:

\begin{bmatrix}e^x&e^{-x}&Cosh\left(x\right)\\e^x&-e^{-x}&Sinh\left(x\right)\\e^x&e^{-x}&Cosh\left(x\right)\end{bmatrix}=0

pues la primera y la última fila del determinante son iguales,y por la propiedades de los determinantes, el resultado es cero.

separador2

Caso de coeficientes constantes: ecuación característica

Cuando los  coeficientes A(x), B(x), C(x) de la ecuación homogénea son constantes A, B, C, que no dependen de x, sabemos que las funciones del tipo y=e^{rx} son soluciones particulares, pues:

y=e^{rx}\Rightarrow Ay''+By'+Cy=Ar^2e^{rx}+Bre^{rx}+Ce^{rx}=e^{rx}\left(Ar^2+Br+C\right)

así que la ecuación homogénea equivale a:

e^{rx}\left(Ar^2+Br+C\right)=0\Rightarrow Ar^2+Br+C=0,

que se reduce a una ecuación algebraica de segundo grado en r, la cual se denomina ecuación característica de la ecuación diferencial lineal homogénea. Resolviéndola, obtenemos dos valores r_{1},r_{2}, y la solución general, por la proposición 1,  será la combinación lineal y=c_{1}e^{r_{1}x}+c_{2}e^{r_{2}x}.

Ejemplo 3: Resolver y''-6y'=0

La ecuación característica se obtiene directamente sustituyendo la "y" por "r" y el orden de la derivada por el exponente correspondiente:  r^2-6r=0; la solución es r^2-6r=0\Rightarrow r_1=0,\;r_2=6, y la solución general es y=c_1e^{0x}+c_2e^{6x}=c_1+c_2e^{6x}.

separador2

Caso de que la ecuación característica tenga solución única

Si la ecuación característica tiene una raíz doble r_{1}=r_{2}=r entonces  no tenemos dos soluciones independientes y_{1},y_{2} ; en tal caso la solución general viene dada por y=\left(c_{1}+c_{2}x\right)e^{rx}.

Ejemplo 4: Resolver y''-2y'+y=0.

La ecuación característica r^{2}-2r+1=0 tiene una única solución r=1, y por tanto la solución general será y=\left(c_{1}+c_{2}x\right)e^{x}.

separador2

Caso de que la ecuación característica tenga soluciones complejas

Cuando la ecuación característica tenga un par de raíces complejas conjugadas r_{1}=a+bi, r_{2}=a-bi, entonces la solució general viene dada por y=e^{ax}\left(c_{1}\cos bx+c_{2}\sin bx\right).

Ejemplo 5: Resolver y''-4y'+5y=0.

La ecuación característica r^{2}-4r+5=0 tiene raíces complejas r_{1}=2+i, r_{2}=2-i, la solución general és y=e^{2x}\left(c_{1}\cos x+c_{2}\sin x\right).

separador2

Ecuaciones lineales de segundo orden no homogèneas

La ecuación lineal no homogènea A(x)y''+B(x)y'+C(x)y=F(x), por ser lineal,  tiene la propiedad de que, si tenemos una solución particular cualquiera y_{p}(x), y si sabemos resolver la ecuación homogénea asociada A(x)y''+B(x)y'+C(x)y=0, siendo y_{h}(x) la solución general
de la homogènea, entonces la solución general de la ecuación no homogénea es y(x)=y_{h}(x)+y_{p}(x). Acabamos de ver que si los coeficientes son constantes, siempre podremos encontrar la solución general de la homogènea.

Ejemplo 6: la función y_{p}(x)=3x es una solución particular de y''+4y=12x. Si resolvemos la ecuación homogènia y''+4y=0 obtenemos  y_{h}(x)=c_{1}\cos x+c_{2}\sin x, así pues la solución general de la no homogénea es y=3x+c_{1}\cos x+c_{2}\sin x.

 separador2

Método de los coeficientes indeterminados

Si en la ecuación lineal completa A(x)y''+B(x)y'+C(x)y=F(x) el término F(x) es un polinomio en x, o bien una exponencial e^{rx}, o bien una función sin(kx) o cos(kx), o bien es una suma de esas funciones, entonces podemos intentar aplicar este método para encontrar una solución particular de la ecuación.

El método consiste en ensayar como solución una combinación lineal parecida a f(x):

  • si f(x) es un polinomio de grado n, ensayamos para y_p un polinomio del mismo grado y coeficientes indeterminados
  • si f(x)=pe^{qx}, ensayamos para y_p una función y_p=Ae^{Bx}
  • si f(x)=sin(kx) , ensayamos para y_p una función y_p=Asin(kx)+Bcos(kx)
  • si f(x) es una suma  de los tipos anteriores, aplicamos la proposición 1 (superposición): encontraremos soluciones particulares para cada tipo por separado, la solución buscada será la suma de las anteriores.

Ejemplo 7: Para encontrar una solución particular de y''+y=2x^3 ensayamos y_p=Ax^3+Bx^2+Cx+D. Derivamos y sustituimos en la ecuación:

\begin{array}{l}y_p=Ax^3+Bx^2+Cx+D;\\y'_p=3Ax^2+2Bx+C;\\y''_p=6Ax+2B;\\y''+y=2x^3\Leftrightarrow\left(6Ax+2B\right)+\left(Ax^3+Bx^2+Cx+D\right)=2x^3\end{array}

Operando e igualando términos:

\begin{array}{l}Ax^3+Bx^2+\left(C+6A\right)x+\left(2B+D\right)=2x^3\Leftrightarrow\\A=2,B=0,\\6A+C=0\Leftrightarrow12+C=0\Leftrightarrow C=-12,\\2B+D=0\Leftrightarrow D=0\end{array}

Hemos obtenido una solución particular: y_p=2x^3-12x.

Ejemplo 8:  Para encontrar una solución particular de y''+y=2x^3-2e^x primero ensayamos y_p=Ax^3+Bx^2+Cx+D y después ensayamos y_p=me^x. La primera solución particular la hemos obtenido en el ejemplo anterior, para la segunda, derivamos y sustituimos:

\begin{array}{l}y_p=me^x\Rightarrow y'_p=y''_p=y;\\y_p''+y_p=-2e^x\Leftrightarrow2y_p=2me^x=-2e^x\Rightarrow m=-1\end{array}

entonces la solución particular será la suma de las obtenidas: y_p=2x^3-12x-e^x.

separador2

Caso de que la solución particular a ensayar sea solución de la homogénea asociada

Cuando en la ecuación lineal completa A(x)y''+B(x)y'+C(x)y=F(x) el término F(x) sea una solución de la homogénea asociada A(x)y''+B(x)y'+C(x)y=0 el método anterior fallará; en ese caso ensayamos y_p=x^a·F(x) siendo a el menor valor entero posible tal que y_p no es solución de la homogénea asociada; tal valor dependerá de la ecuación característica asociada a la homogénea: si tiene raíces simples bastarà con tomar a=1, si tiene raíces dobles tomamos a=2.

Ejemplo 9:  Para encontrar la solución general de y''-y=2-2e^x primero encontramos la solución general de la homogénea: r^2-1=0 y tenemos las raíces r=\pm1\Rightarrow y_h=C_1e^x+C_2e^{-x}. Para encontrar la solución general no podemos ensayar y_p=me^x pues es una solución de la homogénea, y_p=me^x\Rightarrow y''-y=me^x-me^x=0, pero podemos provar y_p=mx^a·e^x, tomando a=1 pues la multiplicidad de las raíces de la homogénea es 1:

\begin{array}{l}y_p=mxe^x\Rightarrow y'_p=me^x+mxe^x=\left(1+x\right)me^x;y''_p=me^x+\left(1+x\right)me^x=\left(2+x\right)me^x;\;\\y''-y=\left(2+x\right)me^x-mxe^x=2me^x=-2e^x\Rightarrow m=-1\end{array}

Tenemos pues la solución general y=y_h+y_p=C_1e^x+C_2e^{-x}-xe^x.

separador2

Método de variación de constantes

El mètodo de coeficientes indeterminados suele ser el más simple para obtener soluciones particulares, pero sólo es útil para ecuaciones con términos independientes de las formas a) polinomios en x, b) exponenciales e^{rx}, c) sen(kx) o cos(kx), d) una combinación lineal de los anteriores. En este apartado vemos un método que, en principio, puede aplicarse siempre, claro que puede generar integrales demasiado complicadas, pero no hay una limitación genérica como sucede con el método de coeficientes indeterminados.

Si tenemos la ecuación lineal de segundo orden completa y''+P(x)y'+Q(x)y=F(x), de la cual hemos obtenido la solución y_h=C_1y_1+C_2y_2 de la homogénea asociada y''+P(x)y'+Q(x)y=0, el método de variación de contantes ensaya una solución y=u_1y_1+u_2y_2 en la que usamos las funciones y_1, y_2 de la y_h pero substituimos las contantes C_1, C_2 por funciones desconocidas u_1, u_2 (de ahí el nombre "variación de constantes"). Entonces tenemos dos incógnitas, u_1, u_2, y una única ecuación que deben cumplir, y''+P(x)y'+Q(x)y=F(x), teniendo más incógnitas que ecuaciones la solución queda indeterminada; podemos añadir una segunda ecuación para eliminar la indeterminación, ¿cuál?

 Derivando la y: y'=u_1'y_1+u_1y_1'+u_2'y_2+u_2y_2'=\left(u_1'y_1+u_2'y_2\right)+\left(u_1y_1'+u_2y_2'\right), en este método se impone que se anule el miembro \left(u_1'y_1+u_2'y_2\right), de forma que la expresión de la derivada de y se simplifica: y'=\left(u_1'y_1+u_2'y_2\right)+\left(u_1y_1'+u_2y_2'\right)=u_1y_1'+u_2y_2'. Volviendo a derivar para obtener la segunda derivada: y''=u_1'y_1'+u_1y_1''+u_2'y_2'+u_2y_2''=\left(u_1'y_1'+u_2'y_2'\right)+\left(u_1y_1''+u_2y_2''\right

Pero tanto y_1 como y_2 son soluciones de la ecuación homogénea asociada, y cumplirán:

    \[y''+P(x)y'+Q(x)y=0\Leftrightarrow y_1''=-P(x)y_1'-Q(x)y_1,\;y_2''=-P(x)y_2'-Q(x)y_2\]

Sustituimos las expresiones de y_1'', y_2'' en la expresión de y'':

\begin{array}{l}y''=\left(u_1'y_1'+u_2'y_2'\right)+\left(u_1y_1''+u_2y_2''\right)=\\\left(u_1'y_1'+u_2'y_2'\right)+u_1\left(-P(x)y_1'-Q(x)y_1\right)+u_2\left(-P(x)y_2'-Q(x)y_2\right)=\\\left(u_1'y_1'+u_2'y_2'\right)-\left(u_1y_1'+u_2y_2'\right)P(x)-\left(u_1y_1+u_2y_2\right)Q(x)=\\\left(u_1'y_1'+u_2'y_2'\right)-y'P(x)-yQ(x)\\\end{array}

Substituyendo en la ecuación completa queda reducida a \begin{array}{l}\left[\left(u_1'y_1'+u_2'y_2'\right)-y'P(x)-yQ(x)\right]+P(x)y'+Q(x)y=F(x)\Leftrightarrow\\u_1'y_1'+u_2'y_2'=F(x)\\\\\end{array}.

En definitiva, el método reduce la ecuación de segundo orden lineal  completa a un sistema de dos ecuaciones lineales de orden 1:

\left.\begin{array}{r}u_1'y_1'+u_2'y_2'=F(x)\\u_1'y_1+u_2'y_2=0\end{array}\right\}

Ejemplo 10: Resolver y''+y=\tan\left(x\right).

La ecuación característica de la homogénea asociada es r^2+1=0 con raíces \pmi, y por tanto la solución general de la homogénea es y_h=C_1\sin\left(x\right)+C_2\cos\left(x\right), identificamos y_1=\sin\left(x\right), y_2=\cos\left(x\right), derivamos ambas: y_1'=\cos\left(x\right), y_2'=-\sin\left(x\right), y escribimos el sistema de ecuaciones:

\left.\begin{array}{r}u_1'\cos\left(x\right)-u_2'\sin\left(x\right)=\tan\left(x\right)\\u_1'\sin\left(x\right)+u_2'\cos\left(x\right)=0\end{array}\right\}

Multiplicamos la 1a ecuación por cos(x) y la 2a por sen(x):

\left.\begin{array}{r}u_1'\cos^2\left(x\right)-u_2'\cos\left(x\right)\sin\left(x\right)=\tan\left(x\right)\\u_1'\sin^2\left(x\right)+u_2'\sin\left(x\right)\cos\left(x\right)=0\end{array}\right\}

Sumamos ambas ecuaciones para eliminar una de las incógnitas, obtenemos una única EDO que en este caso es separable:

\begin{array}{l}u_1'\left(\cos^2\left(x\right)+\sin^2\left(x\right)\right)=\tan\left(x\right)\Leftrightarrow u_1'=\tan\left(x\right)\Leftrightarrow\operatorname du=\tan\left(x\right)\operatorname dx\Leftrightarrow\\u_1=\int\tan\left(x\right)\operatorname dx=-\ln\left(\cos\left(x\right)\right)\end{array}

Substituyendo en alguna de las ecuaciones del sistema, por ejemplo en la primera, obtenemos:

u_1'\sin\left(x\right)+u_2'\cos\left(x\right)=0\Leftrightarrow u_2'=-u_1'\tan\left(x\right)=-\tan^2\left(x\right)

Ahora que ya tenemos u_1, u_2 podemos escribir la solución particular:

\begin{array}{l}y_p=u_1y_1+u_2y_2=-\ln\left(\cos\left(x\right)\right)\cdot\sin\left(x\right)-\tan^2\left(x\right)\cdot\cos\left(x\right)=\\-\sin\left(x\right)\left[\ln\left(\cos\left(x\right)\right)+\tan\left(x\right)\right]\end{array}

y la solución general la obtenemos sumando la de la homogénea y la particular:

y=y_h+y_p=C_1\cos\left(x\right)+C_2\sin\left(x\right)-\sin\left(x\right)\left[\ln\left(\cos\left(x\right)\right)+\tan\left(x\right)\right].

NOTA: De hecho el sistema de dos ecuaciones se puede resolver de forma general, por ejemplo aplicando la regla de Cramer:

\left.\begin{array}{r}u_1'y_1'+u_2'y_2'=F(x)\\u_1'y_1+u_2'y_2=0\end{array}\right\}\Rightarrow u_1'=\frac{\begin{vmatrix}F(x)&y_2'\\0&y_2\end{vmatrix}}{\begin{vmatrix}y_1'&y_2'\\y_1&y_2\end{vmatrix}}=\frac{F(x)y_2}{y_1'y_2-y_2'y_1}

que nos proporciona la primera incógnita:

u_1=\int\frac{F(x)y_2}{y_1'y_2-y_2'y_1}\operatorname dx

y para la segunda tenemos:

u_2'=\frac{\begin{vmatrix}y_1'&F(x)\\y_1&0\end{vmatrix}}{\begin{vmatrix}y_1'&y_2'\\y_1&y_2\end{vmatrix}}=\frac{F(x)y_1}{y_1'y_2-y_2'y_1}\Rightarrow u_2=\int\frac{F(x)y_1}{y_1'y_2-y_2'y_1}\operatorname dx,

con lo cual, para el caso de EDO lineales de segundo orden, obtenemos la denominada fórmula de variación de parámetros, aplicándola, no es necesario seguir los pasos del método:

u_1=\int\frac{F(x)y_2}{y_1'y_2-y_2'y_1}\operatorname dx,\;u_2=\int\frac{F(x)y_1}{y_1'y_2-y_2'y_1}\operatorname dx

De paso, observemos que el determinante del denominador es equivalente al Wronskiano de las funciones y_1, y_2, por tanto siempre que estas soluciones de la homogénea sean linealmente independientes el sistema tendrá solución y existirá la solución particular.

De hecho, la fórmula de variación de parámetros puede generalizarse a EDOs lineales de orden n: si llamamos W_i al determinante que resulta de substituir la columna i-ésima del Wronskiano por el término independiente del sistema \begin{bmatrix}F(x)\\0\\\vdots\\0\end{bmatrix}, resulta la fórmula que nos da cada función u_i(x), que es simplemente:

\boxed{u_i=\int\frac{F(x)W_i}W\operatorname dx}

EDOs lineales de coeficientes constantes de orden superior

Todo lo visto hasta ahora para EDOs lineales de segundo orden puede generalizarse para órdenes superiores.

Ejemplo 11: la ecuación lineal homogénea de tercer orden y'''+y''-2y'=0 tiene la ecuación característica r^3+r^2-2r=0 cuyas raíces podemos encontrar fácilmente, son r=0, 1, -2, entonces la solución general de la ecuación es y=C_1+C_2e^x+C_3e^{-2x}.

Para la ecuación lineal completa de orden n a_ny^{(n}+a_{n-1}y^{(n-1}+\dots a_1y^'+a_0y=F(x) tendremos que resolver primero la homogénea asociada encontrando las n raíces de la ecuación característica, obteniendo n soluciones independientes de la homogénea, y después necesitaremos obtener una solución particular independiente de la completa, para expresar la solución general como y=y_h+y_p=\left(A_ny_{h_n}+A_{h_{n-1}}y+\dots A_1y_{h_1}\right)+y_p La independencia de n soluciones puede comprobarse calculando el Wronskiano, que contiene las funciones y sus derivadas hasta el orden n-1:

\begin{bmatrix}y_1&\cdots&y_n\\\vdots&\vdots&\vdots\\y_1^{(n-1}&&y_n^{(n-1}\end{bmatrix}\neq0

 Ejemplo 12: Resolver la ecuación lineal de cuarto orden y^{(4}-y=3cos(x).

La ecuación característica es r^4-1=0 con raíces r^4-1=0\Leftrightarrow r^2=\pm1\Leftrightarrow r=\pm1,\pm\sqrt{-1}=\pm i, como hemos visto en un apartado anterior, las soluciones independientes de la homogénea son e^x, e^{-x}, cos(x), sin(x), y la solución general de la homogénea es y_h=C_1e^x+C_2e^{-x}+C_3\cos\left(x\right)+C_4\sin\left(x\right). Para la solución particular, como el término independiente es F(x)=3\cos(x) el método de coeficientes indeterminados sugiere y_p=A\cos(x)+B\sin(x) pero en este caso esta función forma parte de la solución de la homogénea, con multiplicidad 1, por tanto tomamos y_p=Ax\cos(x)+Bx\sin(x), derivando cuatro veces obtenemos y_p^{(4}=(4A+Bx)\sin(x)+(Ax-4B)\cos(x), sustituimos en la ecuación completa:

\begin{array}{l}\begin{array}{l}y^{(4}-y=3cos(x)\Leftrightarrow\\(4A+Bx)\sin(x)+(Ax-4B)\cos(x)-Ax\cos(x)+Bx\sin(x)=3cos(x)\Leftrightarrow\end{array}\\4A\sin(x)-4Bcos(x)=3cos(x)\Rightarrow A=0,\;B=-\frac34\end{array},

obtenemos la solución general:

y=C_1e^x+C_2e^{-x}+C_3\cos\left(x\right)+C_4\sin\left(x\right)-\frac34x\cdot\sin\left(x\right)

 separador2

Ecuación de Euler

Todo los que hemos visto hasta ahora son ecuaciones lineales con coeficientes constantes, pues las de coeficientes variables son en general o muy difíciles o imposibles de resolver analíticamente. Uno de los pocos casos resolubles es el de la ecuación de Euler:

a_nx^ny^{(n}+\dots+a_1xy'+a_0y=F(x)

o sea es del tipo en el que las potencias de x disminuyen como los órdenes de derivación de y, incluso cuando no coinciden la potencia y el orden en cada término, ya que entonces podemos multiplicar (o dividir según convenga) toda la ecuación por la potencia de x conveniente. Por ejemplo, la ecuación 3x^2y^{'''}+5xy''-y'\;+\frac yx=e^x se convierte en el tipo Euler multiplicándola por x para obtener 3x^3y^{'''}+5x^2y''-xy'\;+y=xe^x.

Para resolver la homogénea asociada se ensaya y_h=x^r, que proporcionará una ecuación en r. Para encontrar una solución particular se realiza el cambio x=e^t que convierte la ecuación de Euler en otra de coeficientes constantes.

Ejemplo 13: Resolver 3x^3y^{'''}+5x^2y''-xy'\;+y=x.

Probamos y_h=x^r en la homogénea; calculamos las derivadas sucesivas y sustituimos en la ecuación homogénea:

\begin{array}{l}y_h=x^r;\;y'=rx^{r-1};\;y''=r\left(r-1\right)x^{r-2};\;y_h^{(3}=r\left(r-1\right)\left(r-2\right)x^{r-3}\Rightarrow\\3x^3y^{'''}+5x^2y''-xy'\;+y=0\Rightarrow\\3x^3r\left(r-1\right)\left(r-2\right)x^{r-3}+5x^2r\left(r-1\right)x^{r-2}-xrx^{r-1}+x^r=0\Rightarrow\\3r\left(r-1\right)\left(r-2\right)+5r\left(r-1\right)-r+1=0\Leftrightarrow\end{array}

Resolvemos la ecuación: una raíz es inmediato ver que es r=1, dividiendo por r-1 obtenemos la ecuación 3r^2-r-1=0 con raíces \frac{1\pm\sqrt{13}}6. La solución de la homogénea es pues:

y_h=C_1x+C_2x^{\frac{1+\sqrt{13}}6}+C_3x^{\frac{1-\sqrt{13}}6x}

Para encontrar una solución particular de la completa, la transformamos con el cambio x=e^t, primero obtenemos las derivadas de y respecto de la variable t, aplicando la igualdad \frac{\operatorname dy}{\operatorname dx}=\frac{\operatorname dy}{\operatorname dt}\frac{\operatorname dt}{\operatorname dx}, comanzando por y':

y'=\frac{\operatorname dy}{\operatorname dx}=\frac{\operatorname dy}{\operatorname dt}\frac{\operatorname dt}{\operatorname dx}=\frac{\operatorname dy}{\operatorname dt}\left(\frac{\operatorname dx}{\operatorname dt}\right)^{-1}=\frac{\operatorname dy}{\operatorname dt}\left(\frac{\operatorname de^t}{\operatorname dt}\right)^{-1}=\frac{\operatorname dy}{\operatorname dt}e^{-t}

y aplicando la misma técnica reiteradamente:

\begin{array}{l}y''=\frac{\operatorname dy'}{\operatorname dx}=\frac d{dx}\left(\frac{\operatorname dy}{\operatorname dt}e^{-t}\right)=\frac d{dt}\left(\frac{\operatorname dy}{\operatorname dt}e^{-t}\right)\frac{\operatorname dt}{\operatorname dx}=\left(\frac{\operatorname d^2y}{\operatorname dt^2}e^{-t}-\frac{\operatorname dy}{\operatorname dt}e^{-t}\right)e^{-t}=e^{-2t}\left(\frac{\operatorname d^2y}{\operatorname dt^2}-\frac{\operatorname dy}{\operatorname dt}\right);\\y'''=\frac{\operatorname dy''}{\operatorname dx}=\frac d{dx}\left[e^{-2t}\left(\frac{\operatorname d^2y}{\operatorname dt^2}-\frac{\operatorname dy}{\operatorname dt}\right)\right]=\frac d{dt}\left[e^{-2t}\left(\frac{\operatorname d^2y}{\operatorname dt^2}-\frac{\operatorname dy}{\operatorname dt}\right)\right]\frac{\operatorname dt}{\operatorname dx}=\\\left[-2e^{-2t}\left(\frac{\operatorname d^2y}{\operatorname dt^2}-\frac{\operatorname dy}{\operatorname dt}\right)+e^{-2t}\left(\frac{\operatorname d^3y}{\operatorname dt^3}-\frac{\operatorname d^2y}{\operatorname dt^2}\right)\right]e^{-t}=\\e^{-3t}\left[\frac{\operatorname d^3y}{\operatorname dt^3}-3\frac{\operatorname d^2y}{\operatorname dt^2}+2\frac{\operatorname dy}{\operatorname dt}\right]\\\end{array}

Sustituimos en la ecuación original:

\begin{array}{l}3x^3y^{'''}+5x^2y''-xy'\;+y=x\Leftrightarrow\\3e^{3t}e^{-3t}\left[\frac{\operatorname d^3y}{\operatorname dt^3}-3\frac{\operatorname d^2y}{\operatorname dt^2}+2\frac{\operatorname dy}{\operatorname dt}\right]+5e^{2t}e^{-2t}\left(\frac{\operatorname d^2y}{\operatorname dt^2}-\frac{\operatorname dy}{\operatorname dt}\right)-e^t\frac{\operatorname dy}{\operatorname dt}e^{-t}+y=e^t\end{array},

operando, nos queda una ecuación lineal con coeficientes constantes:

3\frac{\operatorname d^3y}{\operatorname dt^3}-4\frac{\operatorname d^2y}{\operatorname dt^2}+y=e^t

ensayamos una solución del tipo Ate^t, ya que Ae^t es una solución de la homogénea:

\begin{array}{l}y=Ate^t;\;\frac{\operatorname dy}{\operatorname dt}=\left(1+t\right)Ae^t;\frac{\operatorname d^2y}{\operatorname dt^2}=\left(2+t\right)Ae^t;\frac{\operatorname d^3y}{\operatorname dt^3}=\left(3+t\right)Ae^t;\\3\frac{\operatorname d^3y}{\operatorname dt^3}-4\frac{\operatorname d^2y}{\operatorname dt^2}+y=e^t\Leftrightarrow\left[3\left(3+t\right)-4\left(2+t\right)+t\right]Ae^t=e^t\Leftrightarrow A=1\end{array}

La solución particular la obtenemos deshaciendo el cambio de variable: y_p=te^t=\ln\left(x\right)\cdot x;\; y la solución general será:

y=y_h+y_p=C_1x+C_2x^{\frac{1+\sqrt{13}}6}+C_3x^{\frac{1-\sqrt{13}}6x}+x\ln\left(x\right)\;.

EDOs de primer orden no lineales, cambios de variable

Introducción

En este post vemos algunos métodos de cambio de variable para resolver ecuaciones diferenciales ordinarias no lineales; en general estas ecuaciones son difíciles de resolver, y muchas veces de hecho no existe una solución analítica, esto es, una función y=f(x) explícita, que sea un polinomio, una función trigonométrica, exponencial, etc. o una combinación de las anteriores, por lo que hay que resolver la ecuación con métodos numéricos. Sólo en algunos casos especiales podemos obtener la solución analítica; incluso en estos casos, puede pasar que no podamos dar la el haz de curvas solución en forma explícita y=f(C,x), sino en forma paramétrica \phi\left(x,y,p,C\right)=0.

¿Comentarios, preguntas? visitad https://www.facebook.com/tallermatematic.

separador2

Ecuaciones sin la y en la que se puede despejar y'

Este es el único caso en el que no aparece un cambio de variable, en principio: cuando la ecuación diferencial no contiene la función y,  F(x,y',y'^2, ...,y^n) = 0, y además puede ponerse en forma de un polinomio en y' tal como

A_0+A_1y'+A_2y'^2+...+A_ny'^n=0,

podemos intentar encontrar las raíces del polinomio, y'_1=F_1(x,y), y'_2=F_2(x,y), ...,y'_n=F_n(x,y), y entonces plantear las n ecuaciones lineales y'=F_1(x,y), y'=F_2(x,y), ...,y'=F_n(x,y).

Ejemplo 1: La ecuación 2y'^2-y'+x=0 es un polinomio de segundo grado en y', que puede resolverse fácilmente:

y'=\frac{1\pm\left(1-8x\right)^{1/2}}4

tenemos dos ecuaciones separables:

\begin{array}{l}\frac{\operatorname dy}{\operatorname dx}=\frac{1+\left(1-8x\right)^{1/2}}4\Rightarrow\int\operatorname dy=\int\frac{1+\left(1-8x\right)^{1/2}}4\operatorname dx\Rightarrow\\y=\int\frac14\operatorname dx+\frac14\frac1{-8}\int-8\left(1-8x\right)^{1/2}\operatorname dx=\frac14\left[x-\frac1{12}\left(1-8x\right)^{3/2}\right]+C\end{array}

la otra ecuación es idéntica salvo en un signo, y resulta y=\frac14\left[x+\frac1{12}\left(1-8x\right)^{3/2}\right]+C. Cada solución tiene su propio haz de curvas; en la figura se representa una curva de cada solución, para C=0.

Dos soluciones de 2y'² - y' + x = 0
Dos soluciones de 2y'² - y' + x = 0

separador2

Ecuaciones en la que se puede despejar y

En este caso de la ecuación F(x,y,y') = 0 podemos despejar y, resultando y=F(x,y'). Se puede intentar la sustitución y'=p que convierte la ecuación en y=F(x,p); derivando respecto de x obtenemos y'=p=F_x+F_p·p' donde la notación F_x, F_p denotan las derivadas parciales de F respecto de x, p. Despejando p' obtenemos p'=\frac{p-F_x}{F_p}=\varphi\left(x,p\right) que es la misma ecuación diferencial pero con respecto a p(x), con la derivada p' despejada, que puede ser más fácil de resolver.

Ejemplo 2: La ecuación x^2y'^2-y=0 permite despejar tanto y' como y; si despejamos y, con la sustitución y'=p, obtenemos x^2p^2=y, derivando respecto de x:

\begin{array}{l}2xp^2+2x^2p\cdot p'=y'=p\Rightarrow2px\left(p+xp'\right)=p\Rightarrow2x\left(p+xp'\right)=1\Rightarrow\\p+xp'=\frac1{2x}\Rightarrow p'+\frac1xp=\frac1{2x^2}\end{array}

que es una ecuación diferencial lineal; la lineal homogénea es p'+\frac1xp=0, de variables separables:

\begin{array}{l}p'+\frac1xp=0\Rightarrow\frac{\operatorname dp}{\operatorname dx}=-\frac1xp\Rightarrow\int\frac{\operatorname dp}p=-\int\frac{\operatorname dx}x\Rightarrow\\\ln\left(p\right)=-\ln\left(x\right)+C\Rightarrow p_h=\frac Cx\end{array}.

La notación p_h significa "solución de la ecuación homogénea" Para encontrar una solución particular de la ecuación completa, utilizamos el método de variación de constantes: p=\frac{C\left(x\right)}x; derivamos y sustituimos en la ecuación lineal:

\begin{array}{l}\begin{array}{l}p'=\frac{C'\left(x\right)\cdot x-C\left(x\right)\cdot1}{x^2}=C'\left(x\right)\frac1x-C\left(x\right)\frac1{x^2};\\p'+\frac1xp=\frac1{2x^2}\Rightarrow\left[C'\left(x\right)\frac1x-C\left(x\right)\frac1{x^2}\right]+\frac1xC\left(x\right)\cdot\frac1x=\frac1{2x^2}\Rightarrow\end{array}\\C'\left(x\right)\frac1x=\frac1{2x^2}\Rightarrow\frac{\operatorname dC}{\operatorname dx}=\frac1{2x}\Rightarrow\int\operatorname dC=\int\frac1{2x}\operatorname dx=\frac12\int\frac1x\operatorname dx\Rightarrow\\C\left(x\right)=\frac12\ln\left(x\right).\end{array}

La solución particular es pues p=\frac{C\left(x\right)}x,\;C\left(x\right)=\frac12\ln\left(x\right)\;\Rightarrow p=\frac{\ln\left(x\right)}{2x} y la solución general se obtiene combinando la solución de la homogénea y la solución particular:

p=C\frac1x+\frac{\ln\left(x\right)}{2x}

Esta igualdad relaciona el parámetro p con la variable independiente x; añadiendo la igualdad x^2p^2=y obtenemos la expresión del haz de curvas solución en paramétricas:

\left\{\begin{array}{l}p=C\frac1x+\frac{\ln\left(x\right)}{2x}\\x^2p^2=y\end{array}\right.

Tal como vienen dadas estas ecuaciones, para determinar los puntos (x,y) de las curvas solución, damos valores a la variable independiente x, a continuación hallamos p(x) con la primera ecuación, y por último hallamos y(x,p) con la segunda ecuación:

Algunos puntos (x,y) obtenidos de las ecuaciones paramétricas
Algunos puntos (x,y) obtenidos de las ecuaciones paramétricas

separador2

Ecuaciones en la que se puede despejar x

Si no podemos despejar y pero sí la x, tendremos x=F(y,y'), la sustitución y'=p convierte la ecuación en x=F(y,p), per ahora consideraremos que p no es función de x, sino de y'; evidentemente, como y' es a su vez función de x, p también ha de ser función de x, pero lo será de forma implícita: p(y)=p\left(y\left(x\right)\right).  Derivando  respecto de x toda la expresiónx=F(y,p) con la regla de la cadena y teniendo en cuenta las funciones implícitas y, p, obtenemos:

\begin{array}{l}\frac{\operatorname dx}{\operatorname dx}=\frac{\partial F(y,p)}{\partial y}\cdot\frac{\operatorname dy}{\operatorname dx}+\frac{\partial F(y,p)}{\partial p}\cdot\frac{\operatorname dp}{\operatorname dy}\cdot\frac{\operatorname dy}{\operatorname dx}\Leftrightarrow\\1=F_y\cdot y'+F_p\cdot p'\cdot y'\\\end{array}

donde la notación F_x, F_p denotan las derivadas parciales de F respecto de x, p. Observar como al derivar F(y,p) respecto de x realmente lo hacemos respecto de las variables y, p, pero siendo éstas últimas a su vez funciones de x, y respectivamente,  multiplicamos por sus derivadas respecto a x, y; por último, en el término \frac{\partial F}{\partial p}F(y,p)\cdot\frac{\operatorname dp}{\operatorname dy}\cdot\frac{\operatorname dy}{\operatorname dx} aparece tambien la derivada de y respecto de x. Esto no es más que la regla de la cadena de las derivadas:

\frac{\operatorname d{}}{\operatorname dx}f\left(g\left(h\left(x\right)\right)\right)=\frac{\operatorname df\left(g\left(h\left(x\right)\right)\right)}{\operatorname dx}\cdot\frac{\operatorname dg\left(h\left(x\right)\right)}{\operatorname dx}\cdot\frac{\operatorname dh\left(x\right)}{\operatorname dx}

Despejamos ahora p'  de la expresión derivada :

\begin{array}{l}p'=\frac{\operatorname dp}{\operatorname dy}=\frac{1-F_y\cdot p}{F_p\cdot p}\\\end{array}

que es una ecuación diferencial en p', con p' despejada y en la que falta la p, o sea el primer caso de este post.

Ejemplo 3: Sea la ecuación diferencial ordinaria no lineal y'y=x. Aplicando el cambio y'=p la ecuación se convierte en py=x; derivando respecto x obtenemos:

\frac{\operatorname dp}{\operatorname dy}\frac{\operatorname dy}{\operatorname dx}y+p\frac{\operatorname dy}{\operatorname dx}=\frac{\operatorname dx}{\operatorname dx}\Rightarrow p'py+p^2=1\Rightarrow p'=\frac{1-p^2}{py}

que es separable: \frac{\operatorname dp}{\operatorname dy}=\frac{1-p^2}{py}\Rightarrow\int\frac p{1-p^2}\operatorname dp=\int\frac1y\operatorname dy, la primera integral es del tipo \int\frac{f'\left(x\right)}{f(x)}=\ln\left(f\left(x\right)\right), la segunda es inmediata:

\begin{array}{l}\int\frac p{1-p^2}\operatorname dp=-\frac12\int\frac{-2p}{1-p^2}\operatorname dp=-\frac12\ln\left(1-p^2\right);\\\int\frac1y\operatorname dy=\ln\left(y\right);\\\int\frac p{1-p^2}\operatorname dp=\int\frac1y\operatorname dy\Rightarrow-\frac12\ln\left(1-p^2\right)=\ln\left(y\right)+C\Rightarrow\\\ln\left(y\cdot\left(1-p^2\right)^{1/2}\right)=C\Rightarrow y\cdot\left(1-p^2\right)^{1/2}=C\Rightarrow\boxed{y=\frac C{\sqrt{1-p^2}}}\\\end{array}

que junto con la igualdad x=py nos da el haz de curvas solución en coordenadas paramétricas \phi\left(x,y,p\right).

NOTA: la ecuación de este último ejemplo es de variables separables y por tanto puede resolverse directamente:

y'y=x\Leftrightarrow\int y\operatorname dy=\int x\operatorname dx\Rightarrow\frac{y^2}2=\frac{x^2}2+C\Rightarrow y^2-x^2=C.

Además, en este caso podemos resolver explícitamente las ecuaciones paramétricas, obteniendo el mismo resultado:

\begin{array}{l}\left.\begin{array}{r}y=\frac C{\sqrt{1-p^2}}\\x=py\end{array}\right\}\Rightarrow p=\frac xy\Rightarrow y=\frac C{\sqrt{1-{\displaystyle\frac{x^2}{y^2}}}}=\frac{Cy}{\sqrt{y^2-x^2}}\Rightarrow\\1=\frac C{\sqrt{y^2-x^2}}\Rightarrow\sqrt{y^2-x^2}\;=C.\end{array},

pero en general, esto no será posible, y tendremos que dejar la solución en paramétricas, como veremos en el siguiente ejemplo.

Ejemplo 4: La ecuación 4y=x^2+(y')^2 es del tipo "se puede despejar y"; probemos el cambio y'=p que la transforma en 4y=x^2+p^2. Derivando respecto de x: 4p=2x+2pp', que es del tipo p'=F(x,p) con F homogénea de grado 0 (ver Funciones homogéneas. Aplicación a las ecuaciones diferenciales):

\begin{array}{l}4p=2x+2pp'\Rightarrow p'=\frac{2p-x}p=F(x,p);\\F(\lambda x,\lambda p)=\frac{2\lambda p-\lambda x}{\lambda p}=\frac{2p-x}p=F(x,p)\end{array},

con el cambio de variable u=p/x se convierte en una ecuación de variables separadas:

\begin{array}{l}u=p/x\Rightarrow u'=\frac{p'x-p}{x^2}\Rightarrow p'=\frac{x^2u'+p}x=\frac{x^2u'+ux}x=xu'+u;\\p'=\frac{2p-x}p\Leftrightarrow xu'+u=\frac{2ux-x}{ux}=\frac{2u-1}u\Rightarrow xu'=\frac{2u-1-u^2}u\Rightarrow\\\frac{\operatorname du}{\operatorname dx}=\frac{2u-1-u^2}{xu}\Rightarrow\int\frac u{2u-1-u^2}\operatorname du=\int\frac1x\operatorname dx\\\\\end{array}

La primera integral es del tipo racional:

\frac u{2u-1-u^2}=\frac u{-\left(u-1\right)^2}=-\frac A{u-1}-\frac B{\left(u-1\right)^2}=-\frac{A\left(u-1\right)+B}{\left(u-1\right)^2}

sustituyendo el valor u=1 llegamos a B=1 y luego a A=1, por tanto:

\int\frac u{2u-1-u^2}\operatorname du=-\int\frac1{u-1}\operatorname du-\int\frac1{\left(u-1\right)^2}\operatorname du=-\ln\left(u-1\right)+\frac1{u-1}

Entonces:

\begin{array}{l}\int\frac u{2u-1-u^2}\operatorname du=\int\frac1x\operatorname dx\Leftrightarrow-\ln\left(u-1\right)+\frac1{u-1}=\ln\left(x\right)+C\Rightarrow\\\ln\left(Cx\right)+\ln\left(u-1\right)=\frac1{u-1}\Rightarrow\ln\left(Cx\left(u-1\right)\right)=\frac1{u-1}\Rightarrow\\Cx=\frac{e^\frac1{u-1}}{u-1}.\end{array}

Si deshacemos el cambio u=p/x obtenemos:

\begin{array}{l}Cx=\frac{e^\frac1{u-1}}{u-1}=\frac1{{\displaystyle\frac px}-1}exp\left(\frac1{\frac px-1}\right)=\frac x{p-x}exp\left(\frac x{p-x}\right)\Rightarrow\\\phi\left(x,p\right)=\frac1{p-x}exp\left(\frac x{p-x}\right)=C\end{array}

que junto con la condición y=(x^2+p^2)/4 nos da el haz de curvas solución dependiente del parámetro p. Pero vemos que la función \phi\left(x,p\right) no permite despejar x, por lo que es de difícil manejo; mejor en este caso dejar las ecuaciones paramétricas en función del parámetro u:

\left.\begin{array}{r}x=C\frac{e^\frac1{u-1}}{u-1}\\y=\frac{x^2+\left(xu\right)^2}4=x^2\frac{1+u^2}4\end{array}\right\}

separador2

Ecuaciones de Bernoulli, Ricatti, Lagrange y Clairaut

Ecuación de Bernoulli

Son de la forma y'+X(x)y=F(x)\cdot y^n y pueden reducirse a lineales con el cambio de variable v(x)=y^{1-n} siempre que n no sea ni 0 ni 1; si n = 0 entonces la ecuación es lineal.

Ejemplo 5: La ecuación dy/dx-3y/2x=2x/y és de Bernoulli, con P(x)=3/2x, Q(x)=2x,  n=-1. El cambio v=y^{2} transforma la ecuación:

y=v^{1/2}\Rightarrow y'=\frac1{2v^{1/2}}\frac{\operatorname dv}{\operatorname dx}\Rightarrow\frac1{2v^{1/2}}v'-\frac{3v^{1/2}}{2x}=\frac{2x}{v^{1/2}}

Multipliquemos todo por 2v^{1/2} para obtener v'-\frac3xv=4x que és lineal. Resolvemos esta ecuación usando la fórmula general (ver el post Ecuaciones diferenciales lineales de 1r orden):

y=e^{-\int X\left(x\right)\operatorname dx}\cdot\left[C+\int F(x)\cdot e^{\int X\left(x\right)\operatorname dx}\operatorname dx\right]

en nuestro caso es:

\begin{array}{l}\begin{array}{l}X(x)=-\frac3x,\;F(x)=4x\;\Rightarrow\int X\left(x\right)\operatorname dx=\int-\frac3x\operatorname dx=-3\ln\left(x\right),\\\int F(x)\cdot e^{\int X\left(x\right)\operatorname dx}=\int4x\cdot e^{-3\ln\left(x\right)}=\int4x\cdot x^{-3}=-4x^{-1},\\v(x)=e^{-\int X\left(x\right)\operatorname dx}\cdot\left[C+\int F(x)\cdot e^{\int X\left(x\right)\operatorname dx}\operatorname dx\right]=\end{array}\\e^{3\ln\left(x\right)}\cdot\left[C-4x^{-1}\right]=Cx^3-4x^2\end{array}

deshacemos el cambio v=y^{2} para obtener y(x)=\pm\sqrt{Cx^3-4x^2}=\boxed{\pm x\sqrt{Cx-4}}.

Ejemplo 6:  La ecuación x\left(dy/dx\right)+6y=3xy^{4/3} no és separable, ni lineal, ni homogénea ni exacta, pero dividiéndola por x resulta dy/dx+6y/x=3y^{4/3}, que es una ecuación de Bernoulli con n=4/3, P(x)=1/x,\: Q(x)=3, aplicando el cambio  v=y^{1-4/3}=y^{-1/3} se transformarà en lineal.

Ecuación de Ricatti

Son de la forma dy/dx=P(x)y^2+Q(x)y+R(x); en el caso particular de que conozcamos una de sus soluciones particulares y_p entonces el cambio y=y_p+1/v la convertirá en lineal.

Ejemplo 7: La ecuación y'+y^2=1+x^2 es de Ricatti; sabiendo que una solución particular es y=x, aplicamos el cambio y=x+1/v, derivamos y sustituimos:

\begin{array}{l}y=x+\frac1v\Rightarrow y'=1-\frac1{v^2}v';\\y'+y^2=1+x^2\Leftrightarrow1-\frac1{v^2}v'+\left(x+\frac1v\right)^2=1+x^2\Rightarrow\\1-\frac1{v^2}v'+x^2+\frac1{v^2}+2x\frac1v=1+x^2\Rightarrow \v'-2xv=1\end{array}

esta última ecuación es lineal en v'.

Ecuación de Lagrange

Son un caso particular de las ecuaciones en las que se pueden despejar y, tienen la forma y=xF(y')+G(y'). Se resuelven como hemos visto, con el cambio y'=p que las transforma en lineales en p, siempre que se cumpla p-F(p)\neq0, condición que no se cumple en las ecuaciones de la forma y=xy'+G(y'), denominadas ecuaciones de Clairaut.

Ecuación de Clairaut

Tienen la forma y=xy'+G(y') y su haz de curvas integral es directamente y=Cx+G(x), fácil de comprobar derivando y sustituyendo en la ecuación.

Ecuaciones diferenciales lineales de 1r orden

Ecuación diferencial lineal de primer orden

Es la que tiene la forma

y'+X(x)y=F(x)

donde X(x) y F(x) son funciones cualesquiera que sólo dependen de la variable independiente x.

Ejemplo 1: Las siguientes ecuaciones de primer orden son lineales

  • y'+y·sin(x)=cos(x)
  • 2y'-5ye^x+x=0, pues operando se convierte en y'+y·(-5/2)e^x = -x
  • x^2y'-e·y+3Ln(x)=0, pues podemos expresar la ecuación como y'-\frac e{x^2}\cdot y=-\frac{3Ln(x)}{x^2}.

En este ejemplo podemos darnos cuenta de que las ecuaciones son lineales en y, no en x. Cuando hay términos no lineales en y, la ecuación deja de ser lineal.

Ejemplo 2: Las siguientes ecuaciones de primer orden no son lineales, pues incluyen términos no lineales en la variable y, que no podemos linealizar:

  • y'-\frac1{x^2}\cdot y^2=x
  • yy'-y\;=2x\;
  • xy'+e^y\;=2xy\;

Solución de la ecuación lineal de primer orden

Ecuaciones lineales homogéneas

Primero nos fijamos en el caso especial de que F(x)=0, denominado ecuación diferencial lineal homogénea: y'+X(x)y=0. En este caso, podemos expresarla en la forma de una ecuación de variables separadas:

\begin{array}{l}y'+X(x)y=0\;\Leftrightarrow\frac{\operatorname dy}{\operatorname dx}=-X(x)y\Leftrightarrow\frac{\operatorname dy}y=-X(x)\operatorname dx\Leftrightarrow\\\int\frac{\operatorname dy}y=-\int X(x)\operatorname dx\Leftrightarrow\ln\left(y\right)=-\int X(x)\operatorname dx+C\Leftrightarrow\\\boxed{y=C\cdot exp\left(-\int X(x)\operatorname dx\right)=C\eta\left(x\right)}\end{array}.

Ejemplo 3: La ecuación diferencial lineal homogénea y'+cos(x)y=0 tiene por solución y=C\cdot exp\left(-\int\cos\left(x\right)\operatorname dx\right)=C\cdot exp\left(-\sin\left(x\right)\right).

NOTA: No debemos confundir las ecuaciones lineales homogéneas, y'+X(x)y=0, con las denominadas ecuaciones diferenciales homogéneas y'=F(x,y) ni con las funciones homogéneas F(x,y), que son tales que cumplen la propiedad F\left(\lambda x,\lambda y\right)=\lambda^nF\left(x,y\right). No tienen nada que ver en absoluto.

Soluciones particulares y solución general

Supongamos ahora que tenemos una ecuación lineal y'+X(x)y=F(x) para la que conocemos una solución particular y_p(x). Entonces, si resolvemos la ecuación homogénea asociada y'+X(x)y=0, con solución y_h(x), se cumple que la solución general de la ecuación "completa" será y=Cy_h+y_p. Es fácil comprobar que es solución sustituyendo en la ecuación completa:

\begin{array}{l}y=Cy_h+y_p\Rightarrow y'=Cy'_h+y'_p\Rightarrow\\y'+X(x)y=\left(Cy'_h+y'_p\right)+X(x)\cdot\left(Cy_h+y_p\right)=\\Cy'_h+X(x)\cdot Cy_h+y'_p+X(x)\cdot y_p=\\C\left(y'_h+X(x)\cdot y_h\right)+\left(y'_p+X(x)\cdot y_p\right)=C\cdot0+F(x)=F(x).\end{array}

Por tanto, en las ecuaciones lineales, si tenemos una solución particular, entonces sabemos hallar la solución general.

Ejemplo 4: La ecuación lineal y'-3xy=xe^{2x^2} tiene la solución particular y_p=e^{2x^2}. Obtenemos la solución de la lineal homogénea asociada, y'-3xy=0, que es y_h=C\cdot exp\left(-\int e^{2x^2}\operatorname dx\right)=C\cdot\frac{-1}{4x}e^{2x^2}. Entonces, la solución general es y=e^{2x^2}+C\cdot\frac{-1}{4x}e^{2x^2}.

Método de variación de constantes

La propiedad anterior con la cual obtenemos la solución general, en la práctica exige el conocimiento de una solución particular de la ecuación, que por lo general no tendremos; no obstante, tenemos un método que puede llevarnos a la obtención de una solución particular conociendo la solución de la ecuación homogénea asociada.

El método de variación de las constantes sustituye la constante C de la solución de la homogénea y_h=C·\eta(x) por una función de x, en principio desconocida: C(x). Se trata entonces de ver si y=C(x)·\eta(x) puede ser una solución particular; para ello, la derivamos y sustituimos en la ecuación lineal completa:

\begin{array}{l}y=C(x)\cdot\eta(x)\Rightarrow y'=C'(x)\cdot\eta(x)+C(x)\cdot\eta'(x);\\y'+X(x)y=F(x)\Leftrightarrow\left(C'(x)\cdot\eta(x)+C(x)\cdot\eta'(x)\right)+X(x)\left(C(x)\cdot\eta(x)\right)=F(x)\Leftrightarrow\\C(x)\cdot\left(\eta'(x)+X(x)\cdot\eta(x)\right)+C'(x)\cdot\eta(x)=F(x)\\\end{array},

teniendo en cuenta que \eta'(x)+X(x)\eta(x)=0 (basta tomar C=1 en la solución de la homogénea y_h=C·\eta(x), nos queda:

\begin{array}{l}C(x)\cdot\left(\cancel{\eta'(x)+X(x)\cdot\eta(x)}\right)+C'(x)\cdot\eta(x)=F(x)\Leftrightarrow\\C'(x)\cdot\eta(x)=F(x)\Leftrightarrow\frac{\operatorname dC}{\operatorname dx}\eta(x)=F(x)\Leftrightarrow\operatorname dC=\frac{F(x)}{\eta(x)}\operatorname dx\Leftrightarrow\\C(x)=\int\frac{F(x)}{\eta(x)}\operatorname dx\\\\\end{array}.

Si podemos resolver esta integral, entonces tenemos una solución particular de la ecuación lineal, y_p=C(x)·\eta(x), que es:

\begin{array}{l}y_p=C(x)\cdot\eta(x)=\int\frac{F(x)}{\eta(x)}\operatorname dx\cdot\eta(x)=e^{-\int X\left(x\right)\operatorname dx}\cdot\int\frac{F(x)}{\eta(x)}\operatorname dx\Leftrightarrow\\\\\boxed{y_p=e^{-\int X\left(x\right)\operatorname dx}\cdot\int F(x)\cdot e^{\int X\left(x\right)\operatorname dx}\operatorname dx.}\end{array}

Ahora también podemos dar la expresión general de la solución completa, obtenida por el método de variación de las constantes:

\begin{array}{l}\begin{array}{l}y=Cy_h+y_p=\\C\cdot e^{-\int X\left(x\right)\operatorname dx}+e^{-\int X\left(x\right)\operatorname dx}\cdot\int F(x)\cdot e^{\int X\left(x\right)\operatorname dx}\operatorname dx\Rightarrow\\\boxed{y=e^{-\int X\left(x\right)\operatorname dx}\cdot\left[C+\int F(x)\cdot e^{\int X\left(x\right)\operatorname dx}\operatorname dx\right]}.\end{array}\\\end{array}

Se puede simplificar un poco esta expresión usando la función solución de la ecuación homogénea \eta\left(x\right), queda:

\boxed{y(x)=\eta\left(x\right)\cdot\left[C+\int F(x)\cdot\eta^{-1}\left(x\right)\operatorname dx\right]}

Por tanto, en el caso de las ecuaciones lineales de primer orden, tenemos una fórmula que las reduce al cálculo de dos integrales.  El método de variación de constantes es general: veremos que se puede aplicar a otras ecuaciones diferenciales.

Ejemplo 5: En un circuito eléctrico que contiene una resistencia R y una inductancia L (abreviadamente, circuito RL) la intensidad de corriente generada al aplicar una tensión V viene dada por la ecuación

RI(t)=V(t)-LI'(t)

 que es una ecuación diferencial lineal en la variable I(t), pues:

RI(t)=V(t)-LI'(t)\Leftrightarrow I'(t)+\frac RLI(t)=V(t)

Si la tensión es constante, V(t)=V, y la ecuación es I'(t)+\frac RLI(t)=V. La ecuación homogénea asociada es I'(t)+\frac RLI(t)=0 que es de variables separables:

\begin{array}{l}\frac{\operatorname dI}{\operatorname dt}=-\frac RLI(t)\Leftrightarrow\frac{\operatorname dI}{I(t)}=-\frac RL\operatorname dt\Leftrightarrow\int\frac{\operatorname dI}{I(t)}=-\int\frac RL\operatorname dt\Leftrightarrow\\\ln\left(I(t)\right)=-\frac RLt+C\Leftrightarrow I(t)=C\cdot exp\left(-\frac RLt\right)=C\eta\left(t\right)\end{array}

La solución general viene dada por

\begin{array}{l}I(t)=\eta\left(t\right)\cdot\left[C+\int X(t)\cdot\eta^{-1}\left(t\right)\operatorname dt\right];\\\int X(t)\cdot\eta^{-1}\left(t\right)\operatorname dt=\int\frac VL\cdot e^{\frac RLt}\operatorname dt=\frac VL\cdot\frac LRe^{\frac RLt}=\frac VRe^{\frac RLt};\\I(t)=e^{-\frac RLt}\cdot\left[C+\frac VRe^{\frac RLt}\right]=\boxed{Ce^{-\frac RLt}+\frac VR}.\\\end{array}

En la práctica se asume
que para t=0 ha de ser I(0)=0, por tanto I(0)=C+\frac VR=0\Rightarrow C=-\frac VR y la expresión de la intensidad queda así:

I(t)=-\frac VRe^{-\frac RLt}+\frac VR=\frac VR\left[1-e^{-\frac RLt}\right]

El término -e^{-\frac RLt} vale 1 para t = 0 y se atenúa ràpidamente, pasando a ser prácticamente nulo: es el denominado régimen transitorio del circuito; después, solo queda el término constante I(t) = V/R (Ley de Ohm): es el régimen estacionario:

Intensidad I(t) en un circuito RL
Intensidad I(t) en un circuito RL

 Unicidad de la solución

El teorema siguiente nos asegura que, en ciertas condiciones, la solución de una ecuación diferencial lineal de primer orden es única, o equivalentemente, que no existen soluciones singulares, toda posible solución está incluida en la solución general.

Teorema 1: Unicidad de la solución de la ecuación diferencial lineal de primer orden.

Dada la ecuación y'+X(x)y=F(x), si las funciones X(x) y F(x) son continuas en un entorno abierto I que contiene a x_0, entonces el problema de valor inicial

y'+X(x)y=F(x), y(x_0)=y_0,

tiene una única solución (la que hemos obtenido, dando un valor adecuado a la constante C).

Ecuaciones diferenciales ordinarias de 1r orden

separador2

Ecuaciones separables

Las ecuaciones diferenciales separables son las de la forma dy / dx = H (x, y) donde H (x, y)  se puede expresar bien como el cociente H (x, y) = f (x) / g (y)  o bien como el producto H (x, y) = f (x) \text {·} g (y).

En el primer caso, dy / dx = f (x) / g (y) \Rightarrow g (y) dy = f (x) dx \Rightarrow \int g (y) dy = \int f (x) dx + C, mientras que en el segundo caso dy/dx=f(x) \text{·}g(y) \Rightarrow g^{-1}(y)dy=f(x)dx \Rightarrow \int g^{-1}(y)dy = \int f(x)dx+C.

Ejemplo 1: Resolver la ecuación diferencial x^{2}\frac{dy}{dx}=\frac{x^{2}+1}{3y^{2}+1}.

Tomemos H(x,y)=\frac{\left(x^{2}+1\right)/x^{2}}{3y^{2}+1}=\frac{f(x)}{g(y)} de forma que \int f(x)dx=\int\frac{x^{2}+1}{x^{2}}dx=x-\frac{1}{x}, \int g(y)dy=\int\left(3y^{2}+1\right)dy=y^{3}+y y por tanto la solución vendrá dada de forma implícita por y^{3}+y=x-\frac{1}{x}+C.

Ejemplo 2: encontrar todas las soluciones de la ecuación dy/dx=2x\sqrt{y-1}.

Hagamos H(x,y) = 2x \sqrt{y-1} = f(x) \text{·}g(y), e integremos cada función: \int f(x)dx=\int2xdx=x^{2}; \int g^{-1}(y)dy=\int\frac{1}{\sqrt{y-1}}dy=2\sqrt{y-1}. Las soluciones son 2\sqrt{y-1}=x^{2}+C\Rightarrow y=\left(\frac{x^{2}+C}{2}\right)^{2}+1 . Observemos que la función y (x) = 1   también es solución de la ecuación, pero no está incluida en la solución general y = \left ( \frac {x ^ {2} + C} {2} \right) ^ {2} +1,  es decir que para todo C tendremos y = \left (\frac {x ^ {2} + C} {2} \right) ^ {2} +1 \neq1. Llamamos a las soluciones no incluidas en la solución general soluciones singulares de la ecuación diferencial.

separador2

Soluciones singulares

En el ejemplo 2 hemos visto que podemos tener soluciones que no están incluidas en la solución general, llamadas soluciones singulares. ¿Cuándo sucederá esto y cómo podemos encontrar esas otras soluciones? Recordemos lo visto en el post Ecuaciones diferenciales ordinarias: introducción, apartado "Haz de curvas y ecuaciones diferenciales de primer orden": la solución general de la ecuación diferencial ordinaria y'=F(x,y) representa un haz de curvas y=f(C,x) tal que cada curva de ese haz cumple la ecuación, o sea que para cada punto (x,y) la pendiente de cualquiera de las curvas y=f(C,x) en ese punto cumple y'=F(x,y).

Curva envolvente de un haz de curvas y solución singular

Dos rectas de un haz (azul y rojo) y curva envolvente del haz (verde)
Dos rectas de un haz (azul y rojo) y curva envolvente del haz (verde)

Dado un haz de curvas F(x,y,C)=0 se denomina curva envolvente del haz a la curva f(x,y) tal que es tangente a todas las curvas del haz.

Ejemplo 3: en la imagen se muestran dos rectas del haz de rectas y=(1-10/C)*x+(10-C), concretamente las dos rectas correspondientes a los valores C=3 en rojo y C=6 en azul, y la curva envolvente del haz, que es tangente a la primera recta en un punto cercano a x=1, y=5 y a la segunda recta cerca de x=3.5, y=1.5.

Recordemos que podemos asociar a cualquier haz de curvas F(x,y,C)=0 una ecuación diferencial y'=f(x,y) tal que el haz es la solución general de la ecuación: para cada valor de C, tenemos una curva del haz tal que su derivada verificay'=f(x,y). Pero hemos visto que la curva envolvente es tangente a todas las curvas del haz, o sea que su derivada en cada punto coincide con la derivada de alguna de las curvas del haz; es por esto que la curva envolvente es también solución de la ecuación diferencial del haz.

Ejemplo 4. Consideremos el haz de parábolas 4y=(x + C)^2, todas las curvas son tangentes al eje X en algún punto; en la imagen se representan tres de las curvas del haz, tangentes a y=0 en los puntos x=-3, -2, -1.

Haz de parábolas 4y=(x + C)²
Haz de parábolas 4y=(x + C)²

¿Cuál es la ecuación diferencial del haz de curvas? Derivamos respecto a x la ecuación del haz: 4y'=2(x + C), y eliminamos la constante C usando esta igualdad y la ecuación del haz:

\left.\begin{array}{r}4y'=2(x+C)\Rightarrow C=2y'-x\\4y=(x\;+\;C)^2\end{array}\right\}\Rightarrow4y=(2y'{)^2=4y'^2}\Rightarrow y=y'^2.

Observemos que la recta y=0, por ser tangente a todas las curvas del haz, tiene la misma derivada que las curvas del haz en cada punto del eje X; por tanto, también verifica la ecuación diferencial del haz. En general: si existe una curva envolvente que es tangente a todas las curvas de un haz, entonces esa curva será una solución singular de la ecuación diferencial del haz.

Cálculo de la envolvente de un haz de curvas

Dado un valor cualquiera x=x_0, la envolvente tendrá la misma derivada en ese punto que alguna curva del haz, esto es, habrá una constante C, dependiente de x=x_0, tal que la derivada de la curva del haz F(x_0,y,C)=0 tendrá el mismo valor en ese punto. Prescindimos de x_0: en general para cada x tendremos una C(x). Derivemos respecto a x la expresión del haz, teniendo en cuenta que y también depende de x, usando la regla de la cadena y la derivación implícita:

F(x,y,C)=0\Rightarrow\frac{\operatorname d{F(x,y,C)}}{\operatorname dx}=\frac{\partial F(x,y,C)}{\partial x}+\frac{\partial F(x,y,C)}{\partial y}\frac{\operatorname dy}{\operatorname dx}=\frac{\partial F(x,y,C)}{\partial x}+\frac{\partial F(x,y,C)}{\partial y}y'=0

Despejamos y' para obtener la derivada en cualquier punto de cualquier curva del haz:

y'=-\frac{\partial F(x,y,C)}{\partial x}\left(\frac{\partial F(x,y,C)}{\partial y}\right)^{-1}.

Pasamos ahora a la curva envolvente; en los puntos de contacto (x,y) entre cada curva y la envolventes las derivadas coinciden con las curvas del haz; en cada valor dado x tendremos determinado un C(x), como el valor y(x) también coincide con el del haz, tendremos que en el punto de contacto se cumple F(x,y,C(x)). Derivamos de nuevo respecto a x:

\begin{array}{l}F(x,y,C\left(x\right))=0\Rightarrow\\\frac{\operatorname d{F(x,y,C\left(x\right))}}{\operatorname dx}=\frac{\partial F(x,y,C)}{\partial x}+\frac{\partial F(x,y,C)}{\partial y}\frac{\operatorname dy}{\operatorname dx}+\frac{\partial F(x,y,C)}{\partial C}\frac{\operatorname dC}{\operatorname dx}\\=\frac{\partial F(x,y,C)}{\partial x}+\frac{\partial F(x,y,C)}{\partial y}y'+\frac{\partial F(x,y,C)}{\partial C}C'=0\end{array}.

Despejamos y':

y'=\left(\frac{\partial F(x,y,C)}{\partial y}\right)^{-1}\left[\frac{\partial F(x,y,C)}{\partial x}+\frac{\partial F(x,y,C)}{\partial C}C'\right].

Las dos expresiones para la derivada y' han de ser iguales:

y'=\left(\frac{\partial F(x,y,C)}{\partial y}\right)^{-1}\left[\frac{\partial F(x,y,C)}{\partial x}+\frac{\partial F(x,y,C)}{\partial C}C'\right]=\left(\frac{\partial F(x,y,C)}{\partial y}\right)^{-1}\left[\frac{\partial F(x,y,C)}{\partial x}\right].

Vemos que ha de ser: \frac{\partial F(x,y,C)}{\partial C}C'=0, que tiene dos soluciones posibles, la trivial C'=0\Rightarrow C=cte y la condición \frac{\partial F(x,y,C)}{\partial C}=0. La primera es una solución sin importancia pues nos dice que la constante C no dependerá
de x (haz de curvas degenerado: todas la curvas son iguales), mientras que la segunda es la condición general que ha de cumplir la envolvente.

Ejemplo 5: encontrar la expresión de la curva envolvente del haz de rectas del ejemplo 3, y=(1-10/C)*x+(10-C).

Lo expresamos en la forma F(x,y,C)=y-(1-10/C)\ast x-(10-C)=0 para derivar respecto a  C e igualar a cero:

\frac{\partial{F(x,y,C)}}{\partial C}=\frac{-10}{C^2}x+1=0\Rightarrow C=\sqrt{10x}.

Sustituimos este valor en la expresión del haz para eliminar C y obtener la ecuación de la envolvente:

F(x,y)=y(1-\frac{10}{\sqrt{10x}})x+(10-\sqrt{10x})=y-10-x+2\sqrt{10x}=0.

La envolvente es y=10+x-2\sqrt{10x}.

Ejemplo 6: calcular la envolvente del haz de curvas 4y=(x + C)^2 del ejemplo 4.

\begin{array}{l}F(x,y,C)=4y-(x+C)^2=0;\\\frac{\partial F}{\partial C}=2(x+C)=0\Rightarrow C=-x;\\F(x,y)=4y-(x-x)^2=4y=0\Rightarrow\boxed{y=0}.\end{array}

separador2

Funciones homogéneas. Aplicación a las ecuaciones diferenciales.

Una función F(x,y) es homogénea de grado n si F(\lambda x,\lambda y)=\lambda^nF(x,y) para todo parámetro \lambda\neq0. Por ejemple, la función F(x,y)=xy-x^2 es homogénea de grado 2, pues F(\lambda x,\lambda y)=\lambda x\lambda y- ( \lambda x)^2=\lambda^2F(x,y).

En el caso especial de que F(x,y) es homogénea de grado 0 tenemos F(\lambda x,\lambda y)=F(x,y), y nos da un método de resolución de ecuaciones diferenciales y'=F(x,y) en los que la F sea homogénea de grado 0: con el cambio de variable u=y/x se convierten en ecuaciones de variables separadas. En efecto:

\begin{array}{l}u=y/x\Leftrightarrow y=ux\Rightarrow y'=u'x+u;\\y'=F(x,y)\Leftrightarrow u'x+u=F(x,ux)=F(1,u)\end{array}

donde hemos usado la x en vez de la \lambda como parámetro de la homogénea. Entonces

\begin{array}{l}\frac{du}{dx}x+u=F(u)\Leftrightarrow\frac{du}{dx}x=F(u)-u\Leftrightarrow\frac{du}{F(u)-u}=\frac{dx}x\Leftrightarrow\\\int\frac{du}{F(u)-u}=\int\frac{dx}x=\ln\left(x\right)+C\end{array}.

Ejemplo 7: Resolver y'=\frac y{\sqrt{x^2+y^2}-x}.

F(\lambda x,\lambda y)=\frac{\lambda y}{\sqrt{\lambda^2x^2+\lambda^2y^2}-\lambda x}=\frac{\lambda y}{\lambda\left(\sqrt{x^2+y^2}-x\right)}=\frac y{\sqrt{x^2+y^2}-x} luego es homogénea de grado 0. Hacemos el cambio u=y/x para conseguir que sea de variables separadas:

\begin{array}{l}y'=u'x+u=\frac{ux}{\sqrt{x^2+u^2x^2}-x}=\frac u{\sqrt{1+u^2}-1};\\\frac{du}{dx}x=\frac u{\sqrt{1+u^2}-1}-u=\frac{u\left(\sqrt{1+u^2}+1\right)}{1+u^2-1}-u=\frac{\sqrt{1+u^2}+1}u-u=\\\frac{\sqrt{1+u^2}+1-u^2}u;\\\frac{udu}{\sqrt{1+u^2}+1-u^2}=\frac{dx}x\\\end{array}.

Luego \int\frac{udu}{\sqrt{1+u^2}+1-u^2}=\int\frac{dx}x=\ln\left(x\right)+C. En la primera integral hacemos el cambio de variable 1+u^2=t^2 que la convierte en:

\int\frac{udu}{\sqrt{1+u^2}+1-u^2}=\int\frac{t\operatorname dt}{t+1-\left(t^2-1\right)}=\int\frac{t\operatorname dt}{-t^2+t+2},

que es de tipo racional; descomponemos en fracciones simples:

\begin{array}{l}\int\frac{t\operatorname dt}{-t^2+t+2}=\int\frac{-1/3}{t+1}+\int\frac{-2/3}{t-2}=-\frac13\ln\left(t+1\right)-\frac23\ln\left(t-2\right)\\=-\frac13\ln\left(\left(t+1\right)\left(t-2\right)^2\right)=\ln\left(\left(t+1\right)\left(t-2\right)^2\right)^{-1/3}\end{array}.

Igualamos las dos integrales, y deshacemos los cambios para obtener la solución general en forma implícita:

\begin{array}{l}\ln\left(\left(t+1\right)\left(t-2\right)^2\right)^{-1/3}=\ln\left(x\right)+C\Leftrightarrow\\Cx=\left(\left(t+1\right)\left(t-2\right)^2\right)^{-1/3}=\left(t^3-3t^2+4\right)^{-1/3}\Leftrightarrow\\\frac C{x^3}=\left(u^2+1\right)^{3/2}-3u^2+1\Leftrightarrow\\C=x^3\left[\left(\left(\frac yx\right)^2+1\right)^\frac32-3\left(\frac yx\right)^2+1\right].\end{array}.

Reducción a homogénea

Las funciones racionales de la forma

F(x,y)=\frac{ax+by+c}{a'x+b'y+c'}

sólo son homogéneas de grado 0 cuando c=c'=0, no obstante si pensamos en el numerador y denominador de la función como si fueran dos rectas ax+by+c=0,\;a'x+b'y+c'=0, podemos probar a encontrar su punto de intersección (x_0,y_0), si existe, entonces el cambio de variable X = x - x_0, Y = y - y_0 (que es una traslación de los ejes xy a los ejes XY, tomando como origen de coordenadas el punto (x_0,y_0)), convierte la función en homogénea de grado 0.

Ejemplo 8:  Resolver y'=\frac{y-x-1}{x+y-1}.

Hallamos la intersección de las rectas:

\left.\begin{array}{r}y-x-1=0\\x+y-1=0\end{array}\right\}y=1,\;x=0

hacemos el cambio X=x, Y=y-1 que transforma la ecuación en:

y'=\frac{y-x-1}{x+y-1}\Leftrightarrow Y'=\frac{Y-X}{X+Y}

que contiene una función f(X,Y) homogénea de grado 0; hacemos el cambio u=Y/X para separar las variables:

\begin{array}{l}u'X+u=\frac{u-1}{1+u}\Leftrightarrow\frac{\operatorname du}{\operatorname dX}X=\frac{u-1}{1+u}-u=\frac{-1-u^2}{1+u}=-\frac{1+u^2}{1+u};\\-\int\frac{1+u}{1+u^2}\operatorname du=\int\frac{\operatorname dX}X=\ln\left(X\right)+C\end{array}.

Descomponemos la primera integral:

\begin{array}{l}\int\frac{1+u}{1+u^2}=\int\frac1{1+u^2}+\int\frac u{1+u^2}=\tan^{-1}\left(u\right)+\frac12\int\frac{2u}{1+u^2}=\\\tan^{-1}\left(u\right)+\frac12\ln\left(1+u^2\right)\end{array}.

Igualamos integrales y deshacemos el cambio:

\begin{array}{l}-\tan^{-1}\left(u\right)-\frac12\ln\left(1+u^2\right)=\ln\left(X\right)+C=\ln\left(CX\right)\Leftrightarrow\\-\tan^{-1}\left(\frac YX\right)=\ln\left(CX\left(1+\left(\frac YX\right)^2\right)\right);\\\ln\left(Cx\left(1+\left(\frac{y-1}x\right)^2\right)\right)+\tan^{-1}\left(\frac{y-1}x\right)=0.\\\end{array}

Caso especial: las rectas son paralelas

En el caso especial de que las dos rectas no tengan intersección, entonces es que son paralelas, con vectores directores proporcionales: \frac ba=\frac{b'}{a'}=r; pero entonces podemos hacer

F(x,y)=\frac{ax+by+c}{a'x+b'y+c}=\frac{a'rx+b'ry+c}{a'x+b'y+c}=\frac{r\left(a'x+b'y\right)+c}{a'x+b'y+c},

con el cambio de variable Y=a'x+b'y se convierte en

\begin{array}{l}F(x,Y)=\frac{rY+c}{Y+c'};\;y=\frac{Y-a'x}{b'}\Rightarrow\frac{\operatorname dy}{\operatorname dx}=\frac1{b'}\left(\frac{\operatorname dY}{\operatorname dx}-a'\right);\\\frac{\operatorname dy}{\operatorname dx}=F(x,y)\Leftrightarrow\frac1{b'}\left(\frac{\operatorname dY}{\operatorname dx}-a'\right)=\frac{rY+c}{Y+c'}\Leftrightarrow\\\frac{\operatorname dY}{\operatorname dx}=b'\frac{rY+c}{Y+c'}+a'\end{array}

que es de variables separadas.

Ejemplo 9: y'=\frac{x+y+1}{x+y-1}.

Las rectas x+y+1=0 y x+y-1 no tienen intersección, son paralelas; hacemos el cambio Y=x+y, Y'=1+y' para obtener:

\begin{array}{l}Y'-1=\frac{Y+1}{Y-1};\;\frac{\operatorname dY}{\operatorname dx}=\frac{Y+1}{Y-1}+1=\frac{Y+1+Y-1}{Y-1}=\frac{2Y}{Y-1};\\\int\frac{Y-1}{2Y}\operatorname dY=\int\operatorname dx;\\\frac12\int\left(1-\frac1Y\right)\operatorname dY=x+C;\\Y-\ln\left(Y\right)=2x+C;\\x+y-\ln\left(x+y\right)=2x+C;\\y=x+\ln\left(x+y\right)+C.\end{array}

separador2

Ecuaciones exactas

Las ecuaciones de la forma M (x, y) + N (x, y) \frac {dy} {dx} = 0, o equivalentemente, M (x, y) dx + N (x, y) dy = 0, diremos que son exactas siempre y cuando M y N  sean funciones continuas y
se verifique la condición \frac {\partial M} {\partial y} = \frac {\partial N} {\partial x}.

La solución vendrá dada de forma implícita por F (x, y) = C donde \frac {\partial F} {\partial x} = M, \: \frac {\partial F} {\partial y} = N.

Ejemplo 10: La ecuación y^{3} dx + 3xy^{2} dy = 0  es exacta, pues N = 3xy^{2}, \frac{\partial N} {\partial x} = 3y^{2}, \: M = y^{3}, \: \frac {\partial M} {\partial y} = 3y^{2}.  Como \frac {\partial F} {\partial x} = M, será F (x, y) = \int MDX = xy^{3} + C (y) y también F (x, y) = \int Ndy = xy^{3} + C (x);  igualando las dos expresiones tenemos F (x, y) = xy^{3}, C (x) = C (y) = 0, así que la solución general es xy^{3} = C.

 

Ejemplo 11: Resolver la ecuación diferencial \left (6xy-y^{3} \right) dx + \left (4y + 3x^{2} -3xy^{2} \right) dy = 0.

Como M = 6xy-y^{3},  N = 4y + 3x^{2} -3xy^{2}, \frac {\partial N} {\partial x} = 6x-3y^{2} = \frac {\partial M} {\partial y}, la ecuación es exacta. Entonces

F (x, y) = \int Ndy = 2y^{2} + 3x^{2} y-xy^{3} + C (x),   F (x, y) = \int MDX = 3x^{2} y-xy^{3} + C(y),

y cuando igualamos las dos expresiones obtenemos C (y) = 2y^{2}.

 

Ecuaciones diferenciales ordinarias: introducción

Motivación

Las ecuaciones que se estudian en el Álgebra, tales como x^{3}+7x^{2}+41=0,
expresan relaciones numéricas estáticas, y su solución son números que cumplen la ecuación. Ahora bien, los fenómenos naturales más interesantes implican relaciones dinámicas y se expresan mejor con relaciones entre cantidades variables, esto es, con ecuaciones que no son estáticas, sino que expresan variaciones y relaciones entre magnitudes cambiantes: son las ecuaciones diferenciales.

Recordemos que la derivada de una función dy/dt=f'(t) se puede interpretar como la proporción de cambio de la variable dependiente y respecto a los cambios en la variable dependiente t. Es por ello las ecuaciones que describen cambios utilizan derivadas de funciones.

Definición 1: una ecuación que contiene una función desconocida y una o más de una de sus derivadas se llama ecuación diferencial. Al resolverla, obtenemos funciones y = f (x) que verifican la ecuación. Si la función sólo tiene una variable independiente, la ecuación se llama ecuación diferencial ordinaria. Si la función depende de dos o más variables, las derivadas serán parciales y la ecuación se llama ecuación en derivadas parciales. El orden de una ecuación diferencial es el de la derivada más alta que aparece en la ecuación.

En este post sólo introduciremos las ecuaciones diferenciales ordinarias, dando algunas definiciones, y resolviendo las más inmediatas.

Ejemplo 1: La variación de la temperatura T de un cuerpo con respecto al tiempo es proporcional a la diferencia (T-A) donde A es la temperatura del medio ambiente (ley de enfriamiento de Newton). La ley física se representa por una ecuación diferencial ordinaria de primer orden (EDO de orden 1):

\frac{dT}{dt}=k\left(T-A\right)

Ejemplo 2: La variación con respecto al tiempo de una población P(t) con índices constantes de nacimiento y mortalidad es proporcional al tamaño de la población, y también es una EDO de orden 1:

\frac{dP}{dt}=kP

Ejemplo 3: La ley de Torricelli establece que la variación respecto al tiempo del volumen V de agua en un depósito que se está vaciando es proporcional a la raíz cuadrada de la profundidad y del agua del depósito:

\frac{dV}{dt}=-ky^{1/2}

Ejemplo 4: La distancia x(t) recorrida en el movimiento acelerado de un cuerpo de masa m sometido a una fuerza variable F(t)  viene dada por una EDO de orden 2:

\frac{d^{2}x}{dt}=\frac{F(t)}{m}

Definición 2: una función y = f (x) se llama solución de una ecuación diferencial si la ecuación se cumple cuando se sustituyen y y sus derivadas por f (x) y sus derivadas, respectivamente. Una solución particular de una ecuación diferencial es cualquier solución que se obtiene asignando valores concretos a las constantes en la solución general. En la práctica, las soluciones particulares se obtienen a partir de condiciones iniciales que proporcionan el valor de la variable dependiente o de alguna de sus derivadas para un valor particular de la variable independiente.

Ejemplo 5:  de la ecuación diferencial y''-y=0 son soluciones: a)  y=sin(x),  b) y=4e^{-x}  c) y=Ce^{-x} para cualquier valor C real. La solución b) es una solución particular obtenida a partir de la solución general c). Aunque es menos obvio, también la solución a) es una solución particular obtenida a partir de la solución general (puede verse usando los desarrollos de Taylor de la función exponencial y de la función seno).

Haz de curvas y ecuaciones diferenciales de primer orden

Otra forma de introducir las ecuaciones diferenciales es desde el punto de vista geométrico. Consideremos la gráfica de la función y=Cx^2 para todos los posibles valores reales de C. En la imagen se representan los valores C=1, 2, 4, 8.

Haz de curvas y=Cx²
Haz de curvas y=Cx²

Se denomina haz de curvas planas al conjunto de todas las curvas que son gráficas de una función general y=F(x,C). Cuando damos a C todos los valores posibles, las curvas generadas van llenando una región R del plano; en el caso de la imagen anterior esa región es todo el semiplano superior x\in\left(-\infty,\infty\right),\;y\in\left[0,\infty\right). En general la región R dependerá del haz.

Nos preguntamos ahora: ¿para cada punto (x,y) del plano existirá un único valor C tal que nos define de forma unívoca la función y=F(x,C) que pase por ese punto? Podemos plantearlo en estos términos: fijando x=x_0 la función F pasa a depender sólo de Cy=F(x_0,C)=f(C); ¿para cada valor de y habrá un único valor de C? Será así siempre que la función f sea estrictamente creciente o decreciente en R, y eso sucederá cuando su derivada no se anule: \frac{\operatorname df}{\operatorname dC}\neq0. En este caso es como si tuviéramos otra función C=\psi(x,y) de dos variables que determina C para cada punto del plano.

En el ejemplo de la imagen anterior, fijando x=x_0 cualquiera obtenemos y=f(C)=Cx_0^2, f'(C)=x_0^2, este valor es siempre distinto de cero excepto en el origen de coordenadas, por tanto para el haz y=Cx^2 tenemos unicidad en el sentido de que dado un \left(x,y\right)\neq\left(0,0\right) existe un único valor C=\psi(x,y)=y/x^2 que determina la curva que pasa por ese punto; por el origen (0,0) en cambio pasan todas las curvas del haz.

Si derivamos la ecuación del haz obtenemos y'=2Cx, entonces, sustituyendo el valor anterior C=\psi(x,y) eliminamos la C para obtener: y'=2Cx=2\left(\frac y{x^2}\right)x=\frac{2y}x, que es la ecuación diferencial del haz de curvas, que es una ecuación de primer orden.

Ejemplo 6: Hallar la ecuación diferencial del haz de curvas plano y=Csin(x).

La región R es la unión de los cuadrantes +X+Y  y -X-Y, en la figura se muestran algunas curvas.

Haz de curvas y=C·Sin(x)
Haz de curvas y=C·Sin(x)

Fijando x=x_0 cualquiera obtenemos y=C\sin\left(x_0\right)=f\left(C\right);f'\left(C\right)=\sin\left(x_0\right). Este valor será cero siempre que x_0=k\mathrm\pi,\;\mathrm k=0,1,\dots. En este conjunto de puntos todas las curvas del haz coinciden, y en el resto de puntos tenemos unicidad: una única curva para cada punto, dada por: C=\frac y{\sin\left(x\right)}=\psi(x,y).

Derivando la ecuación del haz: y'=C\cdot\cos\left(x\right). Sustituimos el valor de C:

\left.\begin{array}{r}y'=C\cdot\cos\left(x\right)\\C=\frac y{\sin\left(x\right)}\end{array}\right\}\Rightarrow y'=\frac{y\cos\left(x\right)}{\sin\left(x\right)}=\frac y{\tan\left(x\right)}

Definición 2: El haz de curvas y=f(C,x) es la solución general de la ecuación diferencial ordinaria; si fijamos un valor de C=C_0, obtenemos una única curva del haz, que denominamos solución particular de la ecuación diferencial ordinaria.

Ejemplo 7: el haz de curvas y=Csin(x) es la solución general de la ecuación diferencial y'=\frac y{\tan\left(x\right)}. Por el punto \left(\frac{\mathrm\pi}2,3\right) pasa una única curva, y=3Sin(x), que es una solución particular de la ecuación diferencial.

Existencia y unicidad de la solución de una ecuación diferencial de primer orden

Hemos visto que para obtener la ecuación diferencial de un haz de curvas hay que derivar y eliminar la constante. Planteamos ahora el problema inverso: ¿dada una ecuación diferencial cualquiera, existe "su" haz de curvas  tal como lo hemos definido? El siguiente teorema nos responde para el caso de ecuaciones de primer orden.

Teorema 1: existencia y unicidad. Si tenemos una ecuación diferencial dada en la forma  y'=f(x,y)  tal que la función f es derivable de todos los órdenes (existen todas las derivadas de cualquier orden) en un entorno de (x_0,y_0), entonces existe una única curva y(x) tal que pasa por el punto (x_0,y_0) y cumple la ecuación y'=f(x,y).

Ecuaciones diferenciales ordinarias de primer orden inmediatas

Vemos en este apartado sólo las EDO de 1r orden más simples de resolver, en otros posts veremos los casos más generales.

Ecuaciones del tipo dy/dx=f\left(x\right)

Son integrables directamente, escribiéndolas como dy=f(x)\Rightarrow y=\int f(x)dx+C.

Ejemplo 8: Resolver la ecuación diferencial y'+x^{2}-1=e^{x}.

La ecuación es equivalente a dy/dx=e^{x}-x^{2}+1\Rightarrow y=\int\left(e^{x}-x^{2}+1\right)dx+C=e^{x}-\frac{1}{3}x^{3}+x+C.

Ecuaciones del tipo dy/dx=ky

Si y = f (x) es una función, las ecuaciones del tipo dy / dx = ky tienen por solución al conjunto de funciones y (x) = Ce ^ {kx} donde C es un número real cualquiera. En general una ecuación diferencial tiene infinitas soluciones.

Ejemplo 9: la solución de la ecuación \frac {dP} {dt} = kp  que establece la evolución de una población P(t) con índices constantes de nacimiento y mortalidad es cualquier función de la forma y (x) = Ce^{kx}.

Ejemplo 10: Supongamos que P(t) es la población de una colonia de bacterias en el tiempo t, que la población en t = 0 es de 1000, y que la población se duplica en una hora. Entonces podemos decir que

1000 = P(0)=C
2000 = P(1)=Ce^{k}

por tanto C=1000 y k=\ln2, luego P(t)=1000e^{x\text{·}\ln2}.

La condición P(0) = 1000 se llama condición inicial pues normalmente el valor t = 0 se toma como el estado inicial. Cuando damos una condición inicial, la solución de la ecuación diferencial ya no tendrá en general infinitas soluciones, sino que tendrá sólo una,  o quizá ninguna si las condiciones son incompatibles. Equivalentemente, al dar una condición inicial pasamos, si existe, de la solución general a la solución particular que cumple esa condición. Así, la condición inicial P (0) = - 1000 no tiene ninguna solución del tipo P(t) = Ce^{kt}.

Ejemplo 11: Para resolver la ecuación dy / dx = \frac{x} {\left (x^{2} +9 \right)^{1/2}} con la condición inicial y(4)=2 hacemos y(x)=\int\frac{x}{\left(x^{2}+9\right)^{1/2}}dx=\left(x^{2}+9\right)^{1/2}+C, como  y(4)=\left(16+9\right)^{1/2}+C=5+C ha de ser C=-3.

separador2

Bibliografía

ECUACIONES DIFERENCIALES - Resumen teórico y colección de ejercicios resueltos, y propuestos.

 

Interpolación polinómica

Introducción

A continuación tenemos una tabla con el censo de población de una província de España en el período de tiempo entre 1920 y 1970 del siglo XX:

Censo histórico, província de Álava, España. Fuente: http://www.ine.es/
Censo histórico, província de Álava, España. Fuente: http://www.ine.es/

Nos preguntamos: ¿podemos usar  estos datos para estimar la población en, digamos, 1965? ¿y para hacer una estimación de la población en el año 1910? Si tuviéramos una función población=f(año) que, para los años dados en la tabla, proporcionase exactamente los valores de la población de la tabla, entonces podríamos usar esa función para aproximar los valores desconocidos. En este post estudiamos cómo conseguir una función polinómica que realizará esa función. Para justificar que realmente se puede hacer, tenemos el siguiente teorema:

Teorema de aproximación polinómica de Weierstrass: Si f:\left[a,b\right]\rightarrow\mathbb{R} es continua, entonces dado un \varepsilon>0 cualquiera, existirá un polinomio P(x) tal que \left|f\left(x\right)-P\left(x\right)\right|<\varepsilon para todo x\in\left[a,b\right].

Polinomio interpolador de Newton

Queremos calcular los coeficientes de un polinomio P(x) tal que coincida exactamente con la tabla de valores dados \left(x_0,y_0\right),\;\left(x_1,y_1\right),\dots,\left(x_n,y_n\right). Supongamos para empezar que solo tenemos dos puntos (x,y) tales como \left(x_0,y_0\right),\left(x_1,y_1\right). El polinomio P(x) más simple que pasa por esos dos puntos es la recta dada por los puntos:

\begin{array}{l}\left.\begin{array}{r}\begin{array}{l}y_0=A+Bx_0\\y_1=A+Bx_1\end{array}\end{array}\right\}\Rightarrow y_1-y_0=B\left(x_1-x_0\right)\Rightarrow B=\frac{y_1-y_0}{x_1-x_0};\\A=y_0-Bx_0=y_0-\frac{y_1-y_0}{x_1-x_0}x_0=\frac{y_0\left(x_1-x_0\right)-\left(y_1-y_0\right)x_0}{x_1-x_0}=y_0\end{array},

o sea la recta de ecuación y=y_0+\frac{y_1-y_0}{x_1-x_0}x. Veamos ahora que esta recta puede ponerse en la forma P(x)=a_0+a_1\left(x-x_0\right), en efecto:

\begin{array}{l}y_0=P(x_0)=a_0+a_1\left(x_0-x_0\right)=a_0\Rightarrow y_0=a_0\\y_1=P(x_1)=a_0+a_1\left(x_1-x_0\right)=y_0+a_1\left(x_1-x_0\right)\Rightarrow\;a_1=\frac{y_1-y_0}{x_1-x_0}\;\end{array}.

Podemos generalizar este polinomio al caso de n+1 puntos (se demuestra por el razonamiento por inducción):

P\left(x\right)=a_0+a_1\left(x-x_0\right)+a_2\left(x-x_1\right)\left(x-x_0\right)+\dots+a_n\left(x-x_n\right)\dots\left(x-x_0\right),

pues tenemos que:

\begin{array}{l}P\left(x_0\right)=a_0+a_1\left(x_0-x_0\right)+a_2\left(x_0-x_1\right)\left(x_0-x_0\right)+\dots+a_n\left(x_0-x_n\right)\dots\left(x_0-x_0\right)=a_0;\\P\left(x_1\right)=a_0+a_1\left(x_1-x_0\right)+a_2\left(x_1-x_1\right)\left(x_1-x_0\right)+\dots+a_n\left(x_1-x_n\right)\dots\left(x_1-x_0\right)=a_0+a_1\left(x_1-x_0\right)\\\cdots\\\end{array}.

De las anteriores igualdades se deduce a_0=y_0,\;a_1=\frac{y_1-y_0}{x_1-x_0}; para el resto de coeficientes es más cómodo en vez de plantear las ecuaciones usar el denominado método de diferencias divididas; para ello introducimos la siguiente notación:

  • diferencias divididas de orden cero: f\left[x_i\right]=f\left(x_i\right)
  • diferencias divididas de orden uno: f\left[x_{i-1},x_i\right]=\frac{f\left[x_i\right]-f\left[x_{i-1}\right]}{x_i-x_{i-1}}
  • diferencias divididas de orden dos: f\left[x_{i-2},x_{i-1},x_i\right]=\frac{f\left[x_{i-1},x_i\right]-f\left[x_{i-2},x_{i-1}\right]}{x_i-x_{i-2}}
  • diferencias divididas de orden k: f\left[x_{i-k},\dots,x_i\right]=\frac{f\left[x_{i-k+1},\dots,x_i\right]-f\left[x_{i-k},\dots,x_{i-1}\right]}{x_i-x_{i-k}}

Con esta notación, el polinomio interpolador de Newton es:

P\left(x\right)=f\left[x_0\right]+f\left[x_0,x_1\right]\left(x-x_0\right)+f\left[x_0,x_1,x_2\right]\left(x-x_0\right)\left(x-x_1\right)+\dots+f\left[x_0,x_1,\dots,x_n\right]\left(x-x_0\right)\left(x-x_1\right)\dots\left(x-x_{n-1}\right).

Ejemplo 1: Calcular el polinomio interpolador de Newton para los datos del censo dados en la introducción, para los años 1950-60-70.

Formamos la tabla de diferencias divididas:

x dif orden 0 dif.orden 1 dif.orden 2
1950 49
1960 69 \frac{(69-49)}{(1960-1950)}=2
1970 133 \frac{(133-69)}{(1970-1960)}=6.4 \frac{6.4-2}{(1970-1950)}=\frac{11}{50}

Los coeficientes del polinomio vienen dados por la diagonal de la tabla: 49, 2, 11/50, de forma que el polinomio es de grado 2: P\left(x\right)=49+2\left(x-1950\right)+\frac{11}{50}\left(x-1950\right)\left(x-1960\right).

Ejemplo 2: Calcular el polinomio interpolador de Newton para los datos del censo dados en la introducción, para los años 1940-60-70

Ampliamos la tabla anterior con una nueva fila y una nueva columna, los cálculos que ya teníamos se mantienen

x dif orden 0 dif.orden 1 dif.orden 2 dif. Orden 3
1940 48
1950 49 \frac{(49-48)}{(1950-1940)}=0.1
1960 69 2 \frac{(2-0.1)}{(1960-1940)}=0.095
1970 133 6.4 \frac{11}{50} \frac{({\displaystyle\frac{11}{50}}-0.095)}{(1970-1940)}=\frac1{240}

Los coeficientes son 48, 0.1, 0.095, 1/240, y el polinomio es de grado 3:

P\left(x\right)=48+0.1\left(x-1940\right)+0.095\left(x-1940\right)\left(x-1950\right)+\frac1{240}\left(x-1940\right)\left(x-1950\right)\left(x-1960\right).

En el siguiente gràfico se representan: (a) los valores reales del censo (puntos azules), (b) el polinomio de grado 2 (rojo), (c) el polinomio de grado 3 (amarillo):

Cens_interpolacio

Observad que el polinomio de grado 2 coincide exactamente con 3 valores reales, los de 1950-60-70, mientras que el de grado 3 lo hace con cuatro valores, 1940-50-60-70. En el caso de 1930 prácticamente coincide, pero no exactamente, necesitaríamos el polinomio de grado 4 para ello, y para tener coincidencia exacta en todos los puntos hay que obtener el de grado 5. Esta es la regla general: con n+1 puntos podemos formar un polinomio de grado n que coincidirá con todos los valores dados. Los cálculos pueden automatizarse fácilmente en una hoja de cálculo.

Ejemplo 3: Obtener el polinomio interpolador de grado 5 de los datos del censo, y dar un valor aproximado de la población en el año 1965.

La tabla completa, obtenida con hoja de cálculo, es:

x  orden 0 orden 1 orden 2 orden 3 orden 4 orden 5
1920 36
1930 40 0.4
1940 48 0.8 0.02
1950 49 0.1 -0.035 -1.83E-003
1960 69 2 0.095 4.33E-003 1.54E-004
1970 133 6.4 0.22 4.17E-003 -4.16E-006 -3.16E-006

En la misma hoja de cálculo disponemos los valores anteriores y formamos una tabla de valores (x,y):

x P(x)
1920 36.0
1930 40.0
1940 48.0
1950 49.0
1960 69.0
1965 95.0
1970 133.0

Observad que todos los valores del polinomio de grado 5 son coincidentes con los de la tabla, hemos incluido el valor interpolado del año 1965.

Limitaciones y mejoras del polinomio interpolador de Newton. Splines.

Cuando el número de puntos a interpolar es grande, caso frecuente en las aplicaciones prácticas, los polinomios serán también de un grado elevado. Aunque el teorema de Weierstrass nos  asegura que podemos aproximar por un polinomio cualquier número n de puntos dados, en la práctica debido a los errores de redondeo y de precisión de los ordenadores esto no es así, de hecho este error aumenta progresivamente con n. Se han propuesto diversas soluciones:

Polinomios interpolantes con distribución de puntos variable: la distancia entre los puntos a interpolar se reduce en los extremos de la distribución. Por ejemplo, para interpolar 100 puntos en el intervalo [-10,10] podemos tomar la distribución de puntos de Chebyshev:

x_k=10\cos\left(\frac{\mathrm{πk}}n\right)=\left\{10\cos\left(0\right),10\cos\left(\frac\pi n\right),\dots10\cos\left(\mathrm\pi\right)\right\}

Distancia entre los puntos sucesivos en la interpolación de chebysev
Distancia entre 100 puntos sucesivos en la interpolación de Chebyshev

Interpolación polinómica por intervalos

Otra solución consiste en evitar el uso de polinomios de grado elevado, interpolando por trozos (subintervalos). Por ejemplo, si tenemos 100 puntos a interpolar, podemos tomarlos de tres en tres, y para cada subconjunto generamos un polinomio de grado 2 que pase por esos puntos. El problema que se presenta entonces está en los puntos de unión de cada subintervalo; si por ejemplo lo aplicamos al problema del censo obtenemos tres polinomios con los intervalos de tiempo 1920-30-40, 1940-50-60 y 1960-70 (éste último con sólo dos puntos es una recta); en los puntos de unión de cada polinomio pueden haber cambios bruscos de dirección que implican la no-derivabilidad en esos puntos:

Interpolación a trozos; en los puntos de unión pueden haber cambios bruscos de dirección
Interpolación a trozos; en los puntos de unión pueden haber cambios bruscos de dirección

Siendo la derivabilidad una condición deseable, este procedimiento puede resultar no aconsejable. Para solucionarlo, se pueden emplear la denominada interpolación por splines, un tipo de interpolación polinómica a trozos que asegura la derivabilidad, y es muy empleada en diseño industrial para aproximar matemáticamente formas diversas.

Resolución aproximada de ecuaciones

  1. Introducción
  2. Aproximaciones, errores, algoritmos
  3. Búsqueda de una raíz de una ecuación
  4. Algoritmo de bisección
  5. Algoritmo de Newton
  6. Algoritmo de la secante
  7. Estratégias para la búsqueda de raíces

separador2

 

Introducción

Consideremos el estudio del crecimiento de poblaciones; el crecimiento de la población puede representarse para intervalos breves de tiempo \triangle t suponiendo que en ese intervalo de tiempo la población crece continuamente con un aumento proporcional al número de individuos presentes en el inicio del intervalo de tiempo; llamando N(t) al número de individuos en el tiempo t entonces \frac{\triangle N(t)}{\triangle t}=\lambda N(t), donde \lambda es la denominada tasa de natalidad. Llevando al límite \triangle t\rightarrow0 obtenemos:

\lim_{\triangle t\rightarrow0}\frac{\triangle N(t)}{\triangle t}=\frac{\operatorname{d}N(t)}{\operatorname{d}t}=N'(t)=\lambda N(t).

Si además tenemos en cuenta que podemos tener inmigración, supongamos para simplificar que a una velocidad de llegada constante v_0, entonces tenemos:

v(t)=N'(t)=\lambda N(t)+v_0,

 que es una ecuación diferencial cuya solución es:

N(t)=N_0e^{\lambda t}+\frac{v_0}\lambda\left(e^{\lambda t}-1\right),

donde N_0 es la población inicial en el instante inicial t=0. Supongamos ahora que esta valor inicial es de 10⁶ individuos, que 435.000 inmigran anualmente y que, pasado un año, la población asciende a 1.564.000 individuos. Si con estos datos queremos hallar la tasa de natalidad \lambda tendremos que resolver la ecuación:

1.564.000=10^6\cdot e^{\lambda}+\frac{435.000}\lambda\left(e^{\lambda}-1\right),

que es una ecuación trascendente, esto es, que no tiene solución analítica. Esta situación se da frecuentemente en las aplicaciones prácticas, y ha de resolverse de forma aproximada empleando los métodos numèricos que vamos a ver en este post. Estos métodos usan algoritmos que se implementan en un ordenador.

Aproximaciones, errores, algoritmos

Cuando utilizamos métodos numéricos para aproximar valores, el valor resultado del cálculo diferirá del valor real exacto: cometeremos un error de aproximación. Hay varias formas de expresar este error:

Definición 1: Si x^* es una aproximación de x, el error absoluto cometido está dado por \left|x-x^\ast\right| y el error relativo por \frac{\left|x-x^\ast\right|}{\left|x\right|}, siempre que x\neq 0.

Ejemplo 1: Cada uno de los vagones de un tren de mercancías tiene un peso aproximado en vacío de 20.000 \pm100Kg , y una longitud de 25 metros \pm0.5. Los errores absolutos cometidos son 100Kg y 0.5m, y los relativos 100/20.000 = 1/200 y 0.5/25 = 1/50 respectivamente. Observemos que el error absoluto del peso es, en las unidades dadas, mucho mayor que el error absoluto en la longitud, en cambio con los errores relativos pasa lo contrario.

Definición 2Si x^* es una aproximación de x con t dígitos significativos, entonces t es el mayor entero no negativo tal que 

\frac{\left|x-x^\ast\right|}{\left|x\right|}<5\cdot10^{-t}

o sea que el error relativo es menor que 5\cdot10^{-t}.

Ejemplo 2: Supongamos que queremos aproximar el valor x=20.000 con cuatro dígitos significativos, entonces ha de ser:

\begin{array}{l}\frac{\left|20000-x^\ast\right|}{\left|20000\right|}<5\cdot10^{-4}\Leftrightarrow\left|20000-x^\ast\right|<20000\cdot5\cdot10^{-4}=10\Leftrightarrow\\19990<x^\ast<20010\end{array}.

Definición 3: Cuando aproximamos un número x por su forma de punto flotante con k dígitos, o notación científica estándar, x=0.d_1d_2...d_k·10^n, el error cometido se llama error de redondeo

Ejemplo 3: Si "cortamos" el número \mathrm\pi=3.14159265\dots a k=4 obtenemos  en notación de punto flotante  \mathrm\pi=0.3141\cdot10^1, si "redondeamos" a 4 decimales como a partir del dígito 5 se suma una unidad al siguiente dígito obtenemos \mathrm\pi=0.3142\cdot10^1. En ambos casos el error cometido será el error de redondeo.

Pérdida de precisión en calculadoras y ordenadores digitales

La representación digital de los números decimales no es exacta, además, las operaciones aritméticas que realizan tampoco son exactas. Hay algunos casos concretos en que los errores pueden ser inaceptables, conviene conocerlos para intentar evitarlos. Son:

  1. división de un número con decimales por un número pequeño
  2. multiplicación de un número con decimales por un número grande
  3. sustracción de dos números casi iguales

Cuando en un algoritmo de cálculo se realizan las operaciones anteriores repetidas veces, el error producido se propaga y puede ir aumentando hasta producir un error total inaceptable. Para evitarlo, es conveniente reformular los cálculos para evitar caer en alguna de las situaciones anteriores.

Ejemplo 4: Resolver x^2+10^4x+1=0 con la máxima precisión permitida por la calculadora (por ejemplo, 8 dígitos). Si usamos la fórmula estándar obtenemos:

x=\frac{-10^4\pm\sqrt{10^8-4}}2=-0.5\cdot10^4\pm4999.9999

Con la solución x_1=-0.5\cdot10^4-4999.9999=-9999.9999 no hay problema, pero con la otra solución sí, pues restamos dos números casi iguales: x_1=-0.5\cdot10^4+4999.9999=-0.1\cdot10^3. En efecto, comprobemos las dos soluciones obtenidas sustituyendo en la ecuación:

\begin{array}{l}x_1^2+10^4x_1+1=9999.9999^2-10^4\cdot9999.9999+1=0\\x_2^2+10^4x_2+1=0.001^2-10^4\cdot0.001+1=-8.999999\end{array}

Vemos que la segunda solución no es aceptable. Para remediarlo, tenemos que evitar restar las dos cantidades tan parecidas al obtener x_2, lo podemos conseguir racionalizando:

\begin{array}{l}x_2=\frac{-b+\sqrt{b^2-4ac}}{2a}\cdot\frac{-b-\sqrt{b^2-4ac}}{-b-\sqrt{b^2-4ac}}\\=\frac{b^2-b^2-4ac}{2a\left(-b-\sqrt{b^2-4ac}\right)}=\frac{-2c}{a\left(-b-\sqrt{b^2-4ac}\right)}\end{array},

sustituimos valores en la calculadora:

\frac{-2c}{a\left(-b-\sqrt{b^2-4ac}\right)}=\frac{-2}{-10^4-\sqrt{10^8-4}}=-1.00000001\cdot10^{-4},

volvemos a comprobar este valor en la ecuación: x_2^2+10^4x_2+1=\left(1.00000001\cdot10^{-4}\right)^2-10^4\cdot1.00000001\cdot10^{-4}+1=0.

Algoritmos de cálculo numérico

Definición 4: Un algoritmo de cálculo numérico es una descripción precisa de cómo realizar una sucesión de cálculos en un orden específico, en vista a obtener un resultado deseado de forma aproximada. Cada paso del algoritmo proporciona un valor aproximado, y en cada paso necesitamos realizar unos cálculos. De esta forma obtenemos, paso a paso, una sucesión de aproximaciones al resultado deseado. Cuando esta sucesión converge al resultado diremos que el algoritmo es convergente. Cada ejecución del conjunto de pasos se llama una iteración del algoritmo.

Ejemplo 5: algoritmo de división para calcular el cociente a/b donde a<b y obtener el resultado con una precisión \varepsilon usando las operaciones DIV (división entera) y MOD (resto de la división entera)

  • Variables de entrada: a, b enteros tales que a<b; \varepsilon>0
  • Variables de salida: resultado
    D
  • Paso 0: hacer D=0, n=1
  • Paso 1: q = 10a DIV b
  • Paso 2: r = 10aMODb
  • Paso 3: D = D + q/10^n
  • Paso 4: si r \neq 0 y q/10^n > \varepsilon, entonces hacer a = r, n=n+1, saltar al paso 1; de lo contrario, escribir "Resultado = " D, finalizar.

Si seguimos los pasos de este algoritmo usando los valores a=3, b=7, \varepsilon=0.001 obtenemos la siguiente secuencia de aproximaciones:

  • Paso 0: hacer D=0, n=1
  • Iteración 1
    • Paso 1: q = 10·3 DIV 7 = 4
    • Paso 2: r = 10·3 MOD 7 = 2
    • Paso 3: D = 0 + 4/10^1 = 0.4, esta es la primera aproximación obtenida
    • Paso 4: como r \neq 0 y q/10^1 = 0.4 > 0.001, hacemos a=r=2, n = 1 + 1 = 2 y saltamos de nuevo al paso 1
  • Iteración 2
    • Paso 1: q = 10·2 DIV 7 = 2
    • Paso 2: r = 10·3 MOD 7 = 6
    • Paso 3: D = 0.4 + 2/10^2 = 0.42, esta es la segunda aproximación obtenida
    • Paso 4: como r \neq 0 y q/10^2 = 0.02 > 0.01, hacemos a=r=6, n = 1 + 1 = 3 y saltamos de nuevo al paso 1
  • Iteración 3
    • Paso 1: q = 10·6 DIV 7 = 8
    • Paso 2: r = 10·6 MOD 7 = 4
    • Paso 3: D = 0.42 + 8/10^3 = 0.428, esta es la tercera aproximación obtenida
    • Paso 4: tenemos r \neq 0 pero  4/10^3 = 0.004 < 0.01, por tanto escribir "Resultado = 0.428", finalizar.

Hemos necesitado tres iteraciones del algoritmo para obtener la precisión deseada. Es sencillo implementar este algoritmo usando una hoja de cálculo, con lo cual podemos obtener una precisión elevada, con 10 iteraciones obtenemos 3/7 con 10 decimales de precisión:

Iteración a b q r D n
0 3 7 0 1
1 2 7 4 2 0.4 2
2 6 7 2 6 0.42 3
3 4 7 8 4 0.428 4
4 5 7 5 5 0.4285 5
5 1 7 7 1 0.42857 6
6 3 7 1 3 0.428571 7
7 2 7 4 2 0.4285714 8
8 6 7 2 6 0.42857142 9
9 4 7 8 4 0.428571428 10
10 5 7 5 5 0.4285714285 11

separador2

Búsqueda de una raíz de una ecuación

Dada una función f(x) uno de los problemas básicos del cálculo numérico es el encontrar un valor x_0 tal que f(x_0)=0, al valor x_0 se le denomina una raíz de f(x). En general, dada una ecuación f(x)=c podemos definir g(x)=f(x)-c y hallando una raíz de g(x) resolvemos la ecuación.

Separación de las raíces de una función

Cuando f(x) tenga más de una raíz, será necesario para aplicar los algoritmos el establecer intervalos [a,b] tales que en cada uno de ellos sólo haya una única raíz; a este proceso se le llama separación de las raíces de la función. Para ello será de gran ayuda tener al menos una idea de la forma de la gráfica de la función. También puede ser de gran ayuda el Teorema de Bolzano (ver Funciones continuas, teorema 4) que reproducimos aquí:

Teorema 1 (de Bolzano)Si f:I\rightarrow\mathbb{R} es una función continua definida en un intervalo cualquiera I=(a,b), y f(a) tiene signo distinto a f(b), entonces existe al menos un punto c intermedio entre a y b tal que f(c)=0.

Una especialización de este teorema, especialmente útil para separar raíces, es el siguiente:

Teorema 2, unicidad de la raíz en un intervaloSi f:I\rightarrow\mathbb{R} es una función continua definida en un intervalo cualquiera I=(a,b),  f(a) tiene signo distinto a f(b), y además f'(x) tiene un signo constante en (a,b), entonces existe un único punto c intermedio entre a y b tal que f(c)=0.

Ejemplo 6: separar las raíces de la función y=x^3+x^2+3.

exemple2_equacions

Es una función continua. La gráfica de la función muestra una única raíz cerca del valor x=-2; sabemos que los polinomios tienen como máximo tantas raíces como indica su grado, tres en este caso. ¿Puede ser que existan más raíces que no vemos? La derivada y'=3x^2+2x para x<-2 toma valores positivos, luego en el intervalo \left(-\infty,-2\right) la función es siempre creciente, además, \lim_{x\rightarrow-\infty}f(x)=-\infty.

Por otra parte \lim_{x\rightarrow\infty}f(x)=\infty, y también para x>1 la derivada es siempre positiva, luego la gráfica no puede volver a cortar el eje X. Así pues hay una única raíz, que separamos diciendo que está, por ejemplo en el intervalo [-2.5,-1.5]. En efecto, f(-2.5)=-51/8<0f(-1.5)=15/8>0, luego por el Teorema de Bolzano ha de existir un x_0\in\left[-2.5,2.5\right] tal que f(x_0)=0.

Algoritmo de bisección

Supongamos que hemos separado una raíz de la función continua f(x) en el intervalo (a,b) de forma que
f(a) tiene distinto signo que f(b). Sea el punto intermedio p_1=(a+b)/2; al calcular f(p_1), si resulta ser que f(p_1)\approx0, o más precisamente, dada la precisión \varepsilon resulta que \left|f(p_1)\right|<\varepsilon, entonces tenemos la raíz aproximada x_0=p_1 y hemos terminado. Caso contrario, o bien f(p_1)>0 o bien f(p_1)<0, también tendremos que  o bien f(p_1)0 tiene el mismo signo que f(a) o bien f(p_1)tiene el mismo signo que f(b). En el primer caso, definimos un nuevo intervalo de separación de la raíz, más estrecho, que será (p_1,b), en el segundo caso el nuevo intervalo será (a,p_1). Repetimos entonces el proceso con el nuevo intervalo, obteniendo en cada iteración un intervalo más estrecho en el cual ha de estar la raíz. Describamos el algoritmo:

Algoritmo de bisección

Para encontrar la raíz p de la función continua f(x) situada en el intervalo (a,b) con una tolerancia \varepsilon:

  • Paso 0: n = 1, introducir el número de iteraciones máximo MAXITER, la tolerancia \varepsilon  y el intervalo (a,b)
  • Paso 1: Si n \leqMAXITER entonces seguir al paso siguiente, de lo contrario parar con el mensaje: "número máximo de iteraciones excedido"
  • Paso 2: Calcular p=(a+b)/2
  • Paso 3: Si \left|f\left(p\right)\right|<\varepsilon o (b-a)/2<\varepsilon entonces parar, el resultado es p, mensaje: "raíz encontrada", en otro caso seguir
  • Paso 4: n = n +1; si f(a) tiene el mismo signo que f(p) entonces hacer a = p, si no hacer b = p. Ir al paso 1.

Ejemplo 7: para hallar la raíz de la función  del ejemplo 6, y=x^3+x^2+3, en el intervalo x_0\in\left[-2.5,2.5\right] aplicamos el algoritmo de bisección con una tolerancia \varepsilon=0.001 y fijando el número máximo de iteraciones en 20; lo implementamos con hoja de cálculo:

n a b p f(a) f(b) f(p) (b-a)/2
1 -2.5 2.5 0 -6.375 24.875 3 2.5
2 -2.5 0 -1.25 -6.375 3 2.609375 1.25
3 -2.5 -1.25 -1.875 -6.375 2.609375 -0.076171875 0.625
4 -1.875 -1.25 -1.5625 -0.076171875 2.609375 1.6267089844 0.3125
5 -1.875 -1.5625 -1.71875 -0.076171875 1.6267089844 0.876739502 0.15625
6 -1.875 -1.71875 -1.796875 -0.076171875 0.876739502 0.4270820618 0.078125
7 -1.875 -1.796875 -1.8359375 -0.076171875 0.4270820618 0.1823334694 0.0390625
8 -1.875 -1.8359375 -1.85546875 -0.076171875 0.1823334694 0.0548227429 0.01953125
9 -1.875 -1.85546875 -1.865234375 -0.076171875 0.0548227429 -0.0102362856 0.009765625
10 -1.865234375 -1.85546875 -1.8603515625 -0.0102362856 0.0548227429 0.0224024495 0.0048828125
11 -1.865234375 -1.8603515625 -1.8627929688 -0.0102362856 0.0224024495 0.0061104308 0.0024414063
12 -1.865234375 -1.8627929688 -1.8640136719 -0.0102362856 0.0061104308 -0.0020560847 0.0012207031
13 -1.8640136719 -1.8627929688 -1.8634033203 -0.0020560847 0.0061104308 0.002028883 0.0006103516

Hemos obtenido el valor aproximado de la raíz p=-1.863, para el cual f(p)=0.004726353, el algoritmo se ha parado en la iteración 13 debido a que (b-a)/2= 0.0006103516 < \varepsilon. Observar que damos el valor de p con tres decimales, pues la precisión viene dada por el intervalo (b-a)/2.

El algoritmo de bisección tiene la propiedad de converger siempre a una solución, pero si el intervalo inicial es demasiado grande, puede ser necesario un número elevado de iteraciones. Dicho en otros términos, tenemos asegurada la convergencia pero ésta puede ser muy lenta. Hay además otro posible fallo: podría ser que se terminara con un valor (b-a)/2 muy pequeño pero en cambio f(p) fuera bastante distinto de cero.

Ejemplo 8: El polinomio f(x)=x^7-3x^6-12x^5+12^4-27x^3+34x^2-14x+18 tiene una raíz en (-4,-3) pues f(-2)=782, f(-1)=-117. Aplicando el algoritmo de bisección con \varepsilon=10^{-5} y fijando el número máximo de iteraciones en 50 obtenemos, después de 16 iteraciones, (b-a)/2=1.526E-005, p=-3.183014, f(p)=-0.03 que no nos da una precisión elevada. Continuando hasta 30 iteraciones conseguimos p=-3.183022, f(p)=2.878E-006.

Algoritmo de Newton

Este algoritmo puede ser explicado usando el desarrollo en serie de Taylor de primer orden de la función f(x) de la cual queremos encontrar una raíz. Sea x_0 una aproximación inicial a la raíz x^*. Sabemos que f\left(x\right)\approx f\left(x_0\right)+f'\left(x_0\right)\cdot\left(x-x_0\right). Igualando a cero y operando:

\begin{array}{l}f\left(x\right)=0\Rightarrow f\left(x_0\right)+f'\left(x_0\right)\cdot\left(x-x_0\right)\approx0\Rightarrow x-x_0\approx-f\left(x_0\right)/f'\left(x_0\right)\Rightarrow\\x\approx x_0-\frac{f\left(x_0\right)}{f'\left(x_0\right)}.\end{array}

Tomando el valor x como siguiente aproximación x_1 y repitiendo el proceso obtenemos el algoritmo de Newton:

Algoritmo de Newton

Para encontrar la raíz x^* de la función continua f(x) situada en el intervalo (a,b) con una tolerancia \varepsilon:

  • Paso 0: n = 1, introducir el número de iteraciones máximo MAXITER, la tolerancia \varepsilon  y el valor inicial x_0
  • Paso 1: Si n \leqMAXITER entonces seguir al paso siguiente, de lo contrario parar con el mensaje: "número máximo de iteraciones excedido"
  • Paso 2: Calcular x_n=x_0-\frac{f\left(x_0\right)}{f'\left(x_0\right)}
  • Paso 3: Si \left|x_n-x_0\right|<\varepsilon entonces parar, el resultado es x_n, mensaje: "raíz encontrada", en otro caso seguir
  • Paso 4:hacer x_0 = x_n, n = n + 1;   Ir al paso 1.

Ejemplo 9: El método de Newton aplicado al polinomio f(x)=x^7-3x^6-12x^5+12^4-27x^3+34x^2-14x+18 con x_0=-1.5  y con \varepsilon=10^{-9},  fijando el número máximo de iteraciones en 10 obtiene:

n X_0 x_n f(x_0) f'(x_0) error
1 -3.5 -3.2799249935 -2204.226 10015.796 0.2200750065
2 -3.2799249935 -3.1950472585 -492.453 5801.9121 0.084877735
3 -3.1950472585 -3.1832336086 -53.912 4563.559 0.0118136499
4 -3.1832336086 -3.1830216906 -0.933 4406.065 0.000211918
5 -3.1830216906 -3.1830216234 -0.0003 4403.272 6.722E-008
6 -3.1830216234 -3.1830216234 -2.837E-011 4403.271 0

Comparando con el ejemplo 8, vemos que el algoritmo de Newton obtiene mejores resultados: con solo 6 iteraciones hemos conseguido un valor muy exacto x=-3.1830216234 tal que f(x)\approx-2.8E-011,  no obstante el método tiene también desventajas.

Propiedades del método de Newton

  • No es válido si resulta que f'(x_n)=0 para alguna de las iteraciones
  • Si el valor inicial x_0 está demasiado lejos del valor exacto x^* el método puede ser divergente: no encontrarà la raíz, de hecho, se alejará de ella
  • Cuando es convergente, es más rápido que el de bisección, excepto si f'(x) toma valores cercanos a 1, en cuyo caso será muy lento (necesitará muchas iteraciones).
  • Para aplicarlo, necesitamos la derivada de f(x); en las aplicaciones prácticas, esto puede suponer un problema.

 

Algoritmo de la secante

 El método de Newton puede verse desde el punto de vista geométrico como la aproximación local de la función por su recta tangente en cada punto x_n. por ejemplo, en la imagen siguiente, para encontrar y=f(x)=0 (línia en negro) por Newton usando como aproximación inicial x_0=3 aproximamos la función por su tangente en ese punto (línea amarilla) y hallamos su punto de corte con el eje X, que es la siguiente aproximación x_1=2. En el método de la secante usamos un intervalo inicial (p_0,p_1) (como en el método de bisección) y aproximamos la función por la recta secante en el intervalo: en la imagen el intervalo es (0,3), el corte de la secante con el eje X es precisamente en x=0 donde f(x)=0: el  método encuentra la raíz en un solo paso, ya que hemos tenido la suerte de escoger como extremo del intervalo la solución de la ecuación.

Recta tangente y recta secante a la gráfica de la función y=f(x)
Recta tangente y recta secante a la gráfica de la función y=f(x)

Caso de que la secante tomada en el intervalo(p_0,p_1) no nos dé el resultado, tomamos una nueva recta secante pasando por los puntos (p_1,p).

Algoritmo de la secante

Para encontrar la raíz x^* de la función continua f(x) situada en el intervalo (a,b) con una tolerancia \varepsilon:

  • Paso 0: n = 1, introducir el número de iteraciones máximo MAXITER, la tolerancia \varepsilon  y el intervalo inicial (p_0,p_1)
  • Paso 1: Si n \leqMAXITER entonces seguir al paso siguiente, de lo contrario parar con el mensaje: "número máximo de iteraciones excedido"
  • Paso 2: Calcular q_0=f(p_0),q_1=f(p_1), p=p_1-q_1\frac{p_1-p_0}{q_1-q_0}
  • Paso 3: Si \left|p-1_1\right|<\varepsilon entonces parar, el resultado es p, mensaje: "raíz encontrada", en otro caso seguir
  • Paso 4: hacer p_0 = p_1, q_0=q_1, p_1=p, q_1=f(p), n = n + 1;   Ir al paso 1.

Ejemplo 10: El método de la secante aplicado al polinomio f(x)=x^7-3x^6-12x^5+12^4-27x^3+34x^2-14x+18 con el intervalo inicial (-4,-3)  y con \varepsilon=10^{-9},  fijando el número máximo de iteraciones en 10 obtiene:

n p0 p1 q0=f(p0) q1=f(p1) p f(p) |p-p1|
1 -4.000000 -3.000000 -10966 609 -3.052613 470.9336 0.05261
2 -3.000000 -3.052613 609.00 470.933 -3.232074 -232.340 0.1794
3 -3.052613 -3.232074 470.933 -232.340 -3.172785 44.3865 0.05928
4 -3.232074 -3.172785 -232.340 44.3865 -3.182295 3.19515 0.0095
5 -3.172785 -3.182295 44.3865 3.19515 -3.183033 -0.04947 0.0007
6 -3.182295 -3.183033 3.1951 -0.04946 -3.183022 5.3810E-005 1.12468E-005
7 -3.183033 -3.183022 -0.0494 5.38105E-005 -3.183022 9.02993E-010 1.22203E-008

El valor obtenido, con 8 cifras de precisión, es p=-3.18302162. Comparando con el método de Newton, ha realizado una iteración más para obtener menor precisión, pero a cambio no ha sido necesario calcular la derivada de la función.

De forma parecida al método de Newton, el de la secante puede diverger si el intervalo inicial es demasiado amplio.

Estratégias para la búsqueda de raíces

Para encontrar las raíces de una función el primer paso será separarlas, como se ha explicado en el apartado Separación de raíces. Para cada intervalo aplicaremos alguno de los algoritmos anteriores. Si algún intervalo es demasiado ámplio los métodos de Newton y de la secante pueden diverger, en ese caso podemos mejorar el intervalo realizando algunas iteraciones del método de bisección, para posteriormente aplicar Newton o la secante. Como ejemplo, resolveremos la ecuación planteada al principio del post, en el problema del crecimiento de las poblaciones.

Ejemplo 11: resolver 1.564.000=10^6\cdot e^{\lambda}+\frac{435.000}\lambda\left(e^{\lambda}-1\right).

Lo expresamos como una función f(\lambda)=1564000-10^6e^\lambda-\frac{435000}\lambda\left(e^\lambda-1\right)=0 de la que queremos buscar raíces. La función es continua excepto quizá en \lambda=0, para ver que pasa en ese punto hacemos el límite:

\begin{array}{l}\lim_{\lambda\rightarrow0}\frac{435000}\lambda\left(e^\lambda-1\right)=\frac00=\lim_{\lambda\rightarrow0}\frac{D_\lambda\left[435000\left(e^\lambda-1\right)\right]}{D_\lambda\left[\lambda\right]}\\=\lim_{\lambda\rightarrow0}\frac{435000e^\lambda}1=435000\end{array}

Por tanto la función es continua en todo su dominio. Tenemos que encontrar un intervalo (a,b) tal que haya un cambio de signo en los extremos del intervalo. Para \lambda=0 hemos visto que la función toma el valor +435000, y se ve claramente que para valores suficientemente grandes de \lambda el signo será negativo. Si probamos con \lambda=1 obtenemos -406829, por tanto el intervalo inicial puede ser (0,1). Para no tener que derivar la función usaremos el método de la secante, exigiendo una tolerancia de 10^{-10} y con un máximo de 10 iteraciones; hay que tener en cuenta que para \lambda=0 obtendremos un error de división por cero, por eso cogemos un valor próximo como extremo inferior del intervalo a=10^{-6}, obtenemos:

n p0 p1 q0=f(p0) q1=f(p1) p f(p) |p-p1|
1 0.000001 1.000000 998999 -406829 0.710613 162480 0.2893
2 1.000000 0.710613 -406829 162481 0.793204 17364 0.0825
3 0.710613 0.793204 162481 17364 0.803086 -866 0.0098
4 0.793204 0.803086 17364 -867 0.802616 4.3290 0.0004
5 0.803086 0.802616 -867 4.32906 0.802619 0.00107 0.000002
6 0.802616 0.802619 4.329056 0.00107 0.80261862306 0 5.7838E-010

La convergencia ha sido rápida porque el intervalo inicial es suficientemente bueno.