GDIDP_AMAT_18_Practica6

264 days ago by etlopez18

Ecuaciones diferenciales

1. 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 $-3$ y $3$ y valores de $x$ entre $-3$ y $3$.

t,x = var('t x') cp0=plot_slope_field(t*x, (t,-3,3), (x,-3,3));cp0.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) sol1=plot(sol0(t,0.3),(t,-3,3),thickness=4,color='green') sol2=plot(sol0(t,0),(t,-3,3),thickness=4,color='blue') sol3=plot(sol0(t,-0.5),(t,-3,3),thickness=4,color='red') (cp0+sol1+sol2+sol3).show(ymin = -3, ymax = 3,figsize=5) 
       

                                
                            

                                

Ecuación logística

Vamos a ver un segundo ejemplo con la llamada ecuación logística:  $x'(t)=T_c(1-x/K)$, donde $T_c$ es la tasa de crecimiento de la población (con recursos ilimitados) y $K$ es la capacidad del sistema, el decir, el máximo número de individuos que admite. 

t,x = var('t x') K = 3 # Capacidad del sistema en miles Tc = 0.7 # Tasa de Incremento de población por unidad de tiempo cp1=plot_slope_field(Tc*(1-x/K)*x, (t,0,5), (x,0,K*2));cp1.show(figsize=5) 
       

                                
                            

                                

Vamos a dibujar varias soluciones de la ecuación anterior sobre el campo de pendientes.

c=var('c') sol1(t,c)=-3*exp((c + t)*7/10)/(1-exp((c + t)*7/10)) sol2(t,c)=3*exp((c + t)*7/10)/(1+exp((c + t)*7/10)) grafico=cp1+plot(sol1(t,1),(t,0,5),thickness=4,color='red')+plot(sol2(t,1/2),(t,0,5),thickness=4,color='blue')+plot(3,(t,0,5),thickness=4,color='green')+plot(0,(t,0,5),thickness=4,color='gray') grafico.show(figsize=5) 
       

                                
                            

                                
 
       

2. Ecuaciones autónomas

Vamos a resolver la ecuación diferencial autónoma de primer orden $x'=x(1-x)$ 

f(x)=x*(1-x) 
       

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

t,x = var('t x') cp2=plot_slope_field(f(x), (t,-3,3), (x,-3,3));cp2.show(figsize=5) 
       

                                
                            

                                

Recordemos que para resolverla despejamos:

\[\frac{dx}{x(1-x)}=dt,\]

e integramos ambos miembros:

\[\int\frac{1}{x(1-x)}\,dx= \int \,dt= t+c\]

h(x)=integrate(1/f(x),x).simplify_full();h # Usamos simplify_full para que nos devuelva la solución simplificada 
       
x |--> log(x/(x - 1))
x |--> log(x/(x - 1))

Para despejar la $x$ usamos solve, que recibe una ecuación (dos funciones con el signo == para indicar que es una ecuación) y despeja la $x$. 

solve(h(x)==t+c,x) 
       
[x == e^(c + t)/(e^(c + t) - 1)]
[x == e^(c + t)/(e^(c + t) - 1)]

Guardamos la solución en una variable.

c=var('c') xs(t,c)= e^(c + t)/(e^(c + t) - 1);show(xs) 
       

                                
                            

                                

Vamos ahora a resolver el problema de valor inicial $x'(t)=x(t)(1-x(t))$, $x(0)=3$. Como tenemos la solución de la ecuación diferencial, basta calcular el valor de la constante para que verifique la condición inicial.

solve(xs(0,c)==3,c) 
       
[c == log(3/2)]
[c == log(3/2)]

Dibujamos la función en el campo de pendientes:

sol=plot(xs(t,log(3/2)),(t,-3,3),thickness=4,color='green') (cp2+sol).show(ymin = -3, ymax = 3,figsize=5) 
       

                                
                            

                                

Nótese que lo que dibuja en medio es la asíntota vertical que tiene la función en x=log(2/3) 

 
       

3. Ecuaciones de variables separadas

Vamos a resolver la ecuación de variables separadas \[x'=\frac{x^2}{t}\] Una primera observación es que no está definida para $t=0$, así que consideraremos sólo los $t$ positivos.

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

                                
                            

                                

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

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

De nuevo utilizamos solve para despejar la $x$.

var('c') solve(ix==it+c,x) 
       
[x == -1/(c + log(t))]
[x == -1/(c + log(t))]
xs(t,c)=-1/(c+log(t)); xs(t,c) 
       
-1/(c + log(t))
-1/(c + log(t))

Vamos a calcular la solución determinada por la condición inicial $x(1)=3$.

solve(xs(1,c)==3,c) 
       
[c == (-1/3)]
[c == (-1/3)]

Mostramos la solución:

c1=-1/3; simplify(xs(t,c1)) 
       
-3/(3*log(t) - 1)
-3/(3*log(t) - 1)

Calculamos hasta donde está definida y la pintamos.

solve((3*log(t) - 1)==0,t) 
       
[t == e^(1/3)]
[t == e^(1/3)]
sol1=plot(simplify(xs(t,c1)),(t,0,e^(1/3)),thickness=4,color='red',ymin=0,ymax=2) (cp3+sol1).show(figsize=5) 
       

                                
                            

                                

 

 

 
       

4. Ecuaciones lineales

Resolvamos la ED $x'=tx+t$, con la c.i. $x(0)=1$. A veces, el comando desolve nos hace el cálculo, como en este caso. En otros casos, expresa la solución de forma implícita.

var('t') x=function('x',t) show(desolve(diff(x,t)-t*x-t==0,x)) 
       
__main__:1: DeprecationWarning: Calling function('f',x) is deprecated.
Use function('f')(x) instead.
See http://trac.sagemath.org/17447 for details.
__main__:1: DeprecationWarning: Calling function('f',x) is deprecated. Use function('f')(x) instead.
See http://trac.sagemath.org/17447 for details.
xs2(t,c)=(c-exp(-1/2*t^2))*exp(1/2*t^2) solve([xs2(0,c)==1],c) 
       
[c == 2]
[c == 2]

La solución es entonces:

xs(t)=xs2(t,2); show(xs) 
       

                                
                            

                                

Para dos condiciones iniciales distintas, tenemos dos curvas distintas:

Curva roja: para la c.i. (0,1)

Curva verde: para la c.i. (0,0), en que la c resultante es 1.

solve([xs2(0,c)==0],c) 
       
[c == 1]
[c == 1]
t,x = var('t x') cp4=plot_slope_field(t*x+t, (t,-2,2), (x,-2,2)); sol1=plot(xs(t),(t,-2,2),thickness=2,color='red') sol2=plot(xs2(t,1),(t,-2,2),thickness=2,color='green') (cp4+sol1+sol2).show(ymin = -2, ymax = 2,figsize=5) 
       

Otra opción es utilizar la fórmula:

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') 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') (cp4+sola+solb).show(ymin = -2, ymax = 2,figsize=5) 
       

                                
                            

                                
 
       

5. 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 aproximar 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.

5.1 Método de Euler

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

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

1 paso

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

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

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 
       
h=1/5.n() 
       
xa 
       
(1.0, 0.0, 0.0, 0.0, 0.0)
(1.0, 0.0, 0.0, 0.0, 0.0)
xa[1]=xa[0]+h*f(0,xa[0]);xa[1] 
       
1.0
1.0
xa[2]=xa[1]+h*f(0+h,xa[1]);xa[2] 
       
1.08
1.08
xa[3]=xa[2]+h*f(0+2*h,xa[2]);xa[3] 
       
1.2464000000000002
1.2464000000000002
xa[4]=xa[3]+h*f(0+3*h,xa[3]);xa[4] 
       
1.5159680000000002
1.5159680000000002
xa[5]=xa[4]+h*f(0+4*h,xa[4]);xa[5] 
       
1.9185228800000003
1.9185228800000003
xa 
       
(1.0, 1.0, 1.08, 1.2464000000000002, 1.5159680000000002,
1.9185228800000003)
(1.0, 1.0, 1.08, 1.2464000000000002, 1.5159680000000002, 1.9185228800000003)

Representación gráfica

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

t,x = var('t x') cp5=plot_slope_field(f(t,x), (t,0,1), (x,0,3),figsize=5); #la solución para la c.i. (0,1) era solucion(t)=2*e^(t^2/2)-1 cp5+line(zip(vector(RDF,range(0,6))*0.2,xa))+plot(solucion(t),(t,0,1),color='red') 
       

5.2 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.

x0=1.0 
       
eulermod(t,x)=x+h/2*(f(t,x)+f(t+h,x+h*f(t,x))) 
       
x1=eulermod(0,x0); x1.n() 
       
1.04000000000000
1.04000000000000
x2=eulermod(0+h,x1); x2.n() 
       
1.16566400000000
1.16566400000000
x3=eulermod(0+2*h,x2); x3.n() 
       
1.39262558720000
1.39262558720000
x4=eulermod(0+3*h,x3); x4.n() 
       
1.75056237504512
1.75056237504512
x5=eulermod(0+4*h,x4); x5.n() 
       
2.28967260055396
2.28967260055396
xa2=[x0,x1,x2,x3,x4,x5] xa2 
       
[1.00000000000000,
 1.04000000000000,
 1.16566400000000,
 1.39262558720000,
 1.75056237504512,
 2.28967260055396]
[1.00000000000000,
 1.04000000000000,
 1.16566400000000,
 1.39262558720000,
 1.75056237504512,
 2.28967260055396]

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),color='purple') cp5+line(zip(vector(RDF,range(0,6))*0.2,xa))+plot(solucion(t),(t,0,1),color='red')+saprox2 
       

5.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 
       
[[0, 1], [0.2, 1.04040266667], [0.4, 1.16657398536], [0.6,
1.39443401578], [0.8, 1.75425283056], [1.0, 2.29743335339]]
[[0, 1], [0.2, 1.04040266667], [0.4, 1.16657398536], [0.6, 1.39443401578], [0.8, 1.75425283056], [1.0, 2.29743335339]]

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=4,color='yellow') (cp5+saprox0+saprox1+saprox2+saprox4+sol1).show(ymin = 0, ymax = 3) 
       

                                
                            

                                

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')