GT_MM_1415_Tema_4_Ecuaciones_Diferenciales

1092 days ago by trinidad1415

Ecuaciones diferenciales

Campo de pendientes

Para pintar el campo de pendientes de una ecuación $x'(t)=f(t,x)$ usamos el comando

plot_slope_field(f(t,x),(t,tmin,tmax),(x,xmin,xmax))

Los valores tmin,tmax,xmin,xmax son respectivamente el menor y mayor valor de t y de x que se dibujarán.

Veamos un ejemplo con $x'(t)=t\ x(t)$ para valores de $t$ entre $-2$ y $2$ y valores de $x$ entre $-2$ y $2$.

t,x = var('t x') sf0=plot_slope_field(t*x, (t,-2,2), (x,-2,2));sf0.show(figsize=5) 
       

Las soluciones de la ecuación anterior son de la forma $x(t)=c e^{\frac{t^2}{2}}$, para distintos valores de la constante $c$. Vamos a dibujar varias soluciones sobre el campo de pendientes.

c=var('c') sol0(t,c)=c*exp(t^2/2) sf0 sol1=plot(sol0(t,0.3),(t,-2,2),thickness=4,color='green') sol2=plot(sol0(t,0),(t,-2,2),thickness=4,color='blue') sol3=plot(sol0(t,-0.5),(t,-2,2),thickness=4,color='red') (sf0+sol1+sol2+sol3).show(ymin = -2, ymax = 2,figsize=5) 
       

Ecuaciones de variables separadas

Vamos a resolver la ecuación de variables separadas \[x'(t)=t\frac{1-x^2}{x}\]

En primer lugar, dibujamos el campo de pendientes para ver cómo serán las soluciones.

f(x)=(1-x^2)/x;g(t)=t; t,x = var('t x') cp2=plot_slope_field(f(x)*g(t), (t,-2,2), (x,-2,2));cp2.show(figsize=5) 
       

$f(x)=(1-x^2)/x$ no está definido cuando $x=0$, así que marcamos esa línea.

Por otra parte, $f(x)=0$ si $x=\pm 1$, así que $1$ y $-1$ son soluciones constantes, que también dibujamos.

soluciones_constantes=plot(-1,(x,-2,2),color='green',thickness=2)+plot(1,(x,-2,2),color='green',thickness=2) no_definido=plot(0,(x,-2,2),color='red',linestyle='--',thickness=2) (cp2+no_definido+soluciones_constantes).show(figsize=5) 
       

Para resolver la ecuación se despeja \[\frac{x}{1-x^2}x'=t\] y se integran ambos miembros:

ix(x)=integrate(x/(1-x^2),x).simplify_full() it(t)=integrate(t,t).simplify_full() ix(x),it(t) 
       

Las soluciones se obtendrán despejando $x$ de 

$$\frac{- log(|x^2-1|)}{2}=\frac{t^2}{2}+c,$$

para cada valor de $c$. En particular, las curvas de nivel de 

$$\frac{- log(|x^2-1|)}{2}-\frac{t^2}{2}$$

contienen a las soluciones.

curvas=contour_plot(-1/2*log(abs(x^2 - 1))- t^2/2,(t,-2,2),(x,-2,2), fill=False, cmap='gray') (cp2+curvas+no_definido+soluciones_constantes).show(figsize=7) 
       

Para calcular todas las soluciones, despejamos la $x$.

var('c') solve(ix==it+c,x) 
       

Tenemos dos soluciones, una para $x>0$ y otra para $x<0$. 

x_negativa(t,c)= -sqrt(e^(-t^2 - 2*c) + 1) x_positiva(t,c)= sqrt(e^(-t^2 - 2*c) + 1) 
       

Vamos a calcular la solución que cumple la condición inicial $x(0)=0.5$. Para ello calculamos el valor de $c$ (utilizamos la solución positiva, pues la $x$ tiene un valor positivo).

solve(x_positiva(0,c)==0.5,c) 
       

Nos han salido dos valores, tomamos el primero (si probamos el segundo coincide el resultado).

c1= log(-2/3*I*sqrt(3)); x1(t)=x_positiva(t,c1).full_simplify() x1.show() 
       

Ahora podemos pintarla.

sol_1=plot(x1(t),(t,-2,2), color='blue', thickness=2) (cp2+curvas+no_definido+soluciones_constantes+sol_1).show(figsize=7) 
       

Vamos a calcular la solución que cumple la condición inicial $x(1)=-0.5$. Para ello calculamos el valor de $c$ (utilizamos la solución negativa, pues la $x$ tiene un valor negativo).

solve(x_negativa(1,c)==-0.5,c) 
       
c2= -1/2*I*pi - 1/2*log(3/4) - 1/2 x2(t)=x_negativa(t,c2).full_simplify() x2.show() 
       

Tenemos que está definida siempre que $3e-4e^{t^2}<0$. Calculamos en qué punto la expresión anterior se anula:

solve(3*e-4*exp(t^2)==0,t) 
       

Es decir, que nuestra solución estará definida a partir de $t=\sqrt{\log(3 e/4)}$. La dibujamos.

sol_2=plot(x2(t),(t,sqrt(log(3/4*e)),2), color='purple', thickness=2) (cp2+curvas+no_definido+soluciones_constantes+sol_1+sol_2).show(figsize=7) 
       

Ecuaciones lineales

Si la ecuación lineal es $x'(t)=a(t) x(t)+b(t)$, y definimos

\[A(t)=\int_{t_0}^t a(s)\, ds\] tenemos \[x(t)=x(t_0)e^{A(t)}+e^{A(t)}\int_{t_0}^t b(s) e^{-A(s)}\,ds. \]

Vamos a utilizarlo para calcular las soluciones de la ecuación diferencial $x'=tx+t$.

s=var('s') a(t)=t b(t)=t A(t,t0)=integrate(s,s,t0,t) xs(t,t0,x0)=x0*exp(A(t,t0))+exp(A(t,t0))*integrate(b(s)*exp(-A(s,t0)),s,t0,t);show(simplify(xs)) 
       

Ahora, para resolver cada problema de valor inicial, únicamente tenemos que sustituir $t_0$ y $x_0$ en la función anterior.

t,x = var('t x') cp3=plot_slope_field(t*x+t, (t,-2,2), (x,-2,2)); sola=plot(simplify(xs(t,0,1)),(t,-2,2),thickness=2,color='red') solb=plot(simplify(xs(t,0,0)),(t,-2,2),thickness=2,color='green') (cp3+sola+solb).show(ymin = -2, ymax = 2) 
       

Resolución numérica

Vamos a resolver numéricamente la ecuación diferencial $x'=tx+t$ determinada por la condición inicial $x(0)=1$ e intentar obtener el valor de la solución en $t=1$. Esta ecuación ya la hemos resuelto en el apartado anterior, así que podremos comparar las soluciones obtenidas con la solución exacta.

Método de Euler

Definimos la función $f(t,x)=tx+t$.

f(t,x)=t * x+ t 
       

1 paso

Podemos integrar mediante un paso del método de Euler.

xaprox=1+f(0,1);xaprox 
       

5 pasos

Lo normal será dividir el segmento en varios intervalos y aplicar el método en varios pasos. Vamos a dividir el intervalo $[0,1]$ en cinco pasos y utilizar un vector xa para guardar las soluciones.

xa=vector(RDF,6);xa[0]=1; 
       
xa[1]=xa[0]+0.2*f(0,xa[0]);xa[1] 
       
xa[2]=xa[1]+0.2*f(0.2,xa[1]);xa[2] 
       
xa[3]=xa[2]+0.2*f(0.4,xa[2]);xa[3] 
       
xa[4]=xa[3]+0.2*f(0.6,xa[3]);xa[4] 
       
xa[5]=xa[4]+0.2*f(0.8,xa[4]);xa[5] 
       

Representación gráfica

Vamos a representar el campo de pendientes, la solución aproximada por un paso de Euler (verde), la solución aproximada por cinco pasos de Euler (azul) y la solución exacta (rojo).

t,x = var('t x') cp4=plot_slope_field(f(t,x), (t,0,1), (x,0,4)); saprox0=line(([0,1],[1,xaprox]),thickness=2,color='green') saprox1=line(zip(vector(RDF,range(0,6))*0.2,xa),thickness=2,color='blue') sol1=plot(simplify(xs(t,0,1)),(t,0,1),thickness=2,color='red') (cp4+saprox0+saprox1+sol1).show(ymin = 0, ymax = 3) 
       

Método de Euler modificado

Para aplicar el Método de Euler modificado procedemos análogamente solo cambiando la fórmula para obtener el valor de la aproximación.

xa2=vector(RDF,6) xa2[0]=1; 
       
xa2[1]=xa2[0]+0.1*(f(0,xa2[0])+f(0.2,xa2[0]+0.2*f(0,xa2[0])));xa2[1] 
       
xa2[2]=xa2[1]+0.1*(f(0.2,xa2[1])+f(0.4,xa2[1]+0.2*f(0.2,xa2[1])));xa2[2] 
       
xa2[3]=xa2[2]+0.1*(f(0.4,xa2[2])+f(0.6,xa2[2]+0.2*f(0.4,xa2[2])));xa2[3] 
       
xa2[4]=xa2[3]+0.1*(f(0.6,xa2[3])+f(0.8,xa2[3]+0.2*f(0.6,xa2[3])));xa2[4] 
       
xa2[5]=xa2[4]+0.1*(f(0.8,xa2[4])+f(1.0,xa2[4]+0.2*f(0.8,xa2[4])));xa2[5] 
       

Lo añadimos a la gráfica anterior. El método de Euler modificado es la gráfica en color morado que prácticamente se confunde con la de la solución, en rojo.

saprox2=line(zip(vector(RDF,range(0,6))*0.2,xa2),thickness=2,color='purple') (cp4+saprox0+saprox1+saprox2+sol1).show(ymin = 0, ymax = 3) 
       

Método de Runge-Kutta

Repitiendo los pasos anteriores cambiando la fórmula obtendríamos el método de Runge-Kutta. Sin embargo vamos a aprovechar que dicho método está implementado en Sage. Para ello usamos desolve_rk4, que recibe la función $f(t,x)$ (en nuestro caso $tx+t$), la incógnita ($x$), la condición inicial ($[0,1]$), el punto final y el paso y aplica el método.

x,t=var('x t') sap4=desolve_rk4(t*x+t,x,ics=[0,1],end_points=1,step=0.2);sap4 
       

Dibujamos todo junto. Podemos ver que en el ejemplo el error cometido por el método de Runge-Kutta (en amarillo la solución obtenida por este método) es muy pequeño.

saprox4=line(sap4,thickness=2,color='yellow') (cp4+saprox0+saprox1+saprox2+saprox4+sol1).show(ymin = 0, ymax = 3) 
       

Errores y ecuaciones en dimensión superior

Soluciones a tiempo grande

Los errores de la resolución numérica (por cualquiera de los métodos) son acumulativos, por lo que para tiempos largos los métodos pueden fallar. Vamos a ilustrarlo con el siguiente ejemplo. Tomamos el problema de valor inicial \[x'(t)=x(t)-\sin(t)+\cos(t),\quad x(0)=0,\] cuya solución es $x(t)=\sin(t)$. En la siguiente gráfica se muestra en azul la solución exacta y en amarillo la aproximada. Podemos ver cómo ambas se separan cuando el tiempo aumenta.

x,t=var('x t') sap5=desolve_rk4(x-sin(t)+cos(t),x,ics=[0,0],end_points=14,step=0.2); line(sap5,thickness=4,color='yellow')+plot(sin(t),(t,0,14),thickness=4,color='blue') 
       

Vaguada

En este último ejemplo, calcularemos la vaguada partiendo de un collado. Para ello supondremos que las cotas de cada punto vienen dadas por una función $f(x,y)$.

f(x,y)=-(x^2*y^2+y^4+x^2/2+x*y) niveles=contour_plot(f(x,y),(x,-0.6,0.6),(y,-0.6,0.6), fill=False, cmap='hsv',labels=True, label_fmt="%1.2f", label_colors='black',contours=[0,-0.1,-0.02,-0.05,0.02,0.03,-0.2,-0.3,-0.4,-0.6,-0.8]) niveles.show(figsize=5) 
       

Cargamos el Runge-Kutta 4 para sistemas de ecuaciones.

from sage.calculus.desolvers import desolve_system_rk4 x,y,t=var('x y t') 
       

La ecuación que sigue la vaguada viene dada por el gradiente. Lo calculamos.

g=f.diff(); show(g(x,y)) 
       
J=g.diff();show(J(x,y)) 
       

Buscaremos las direcciones en las que sale del collado. Vienen dadas por los autovectores de la matriz jacobiana.

K=J(0,0) CI1=vector(RDF,K.eigenvectors_right()[0][1][0]) CI2=vector(RDF,K.eigenvectors_right()[1][1][0]) CI1,CI2 
       

Resolvemos las ecuaciones.

P1=desolve_system_rk4(-g(x,y),[x,y],ics=[0,CI1[0]/1000,CI1[1]/1000],ivar=t,end_points=4) Q1=[ [j,k] for i,j,k in P1] LP1=list_plot(Q1,plotjoined=True,color='black',linestyle="--") 
       
P2=desolve_system_rk4(-g(x,y),[x,y],ics=[0,-CI1[0]/1000,-CI1[1]/1000],ivar=t,end_points=4) Q2=[[j,k] for i,j,k in P2] LP2=list_plot(Q2,plotjoined=True,color='black',linestyle="--") 
       

Lo pintamos.

LP1+LP2+niveles