GDIDP_AMAT_1415_EjerciciosyPractica6

1013 days ago by evalsanjuan

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

Ejercicio 1.

a) Dibujar el campo de pendientes de la ecuación diferencial $x'(t)=\cos(t) x(t)$ para $t$ entre $0$ y $2 \pi$ y $x$ entre $-2$ y $2$. 

b) Añadir a dicho campo las gráficas de las funciones $e^{\sin(t)}/2$ y $-1-\sin(t)$. 

c) Discutir cuál de ellas es una posible solución.


 
       

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{x'(t)}{x(t)(1-x(t))}=1,\]

e integramos ambos miembros:

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

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'(t)=\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) 
       

 

 

Ejercicio 2: Resolver la ecuación $x'=2x/t$, pintar el campo de pendientes de la misma, y varias soluciones particulares.

 
       

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)) 
       
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$.

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 
       
xa 
       
(1.0, 0.0, 0.0, 0.0, 0.0, 0.0)
(1.0, 0.0, 0.0, 0.0, 0.0, 0.0)
xa[1]=xa[0]+0.2*f(0,xa[0]);xa[1] 
       
1.0
1.0
xa[2]=xa[1]+0.2*f(0.2,xa[1]);xa[2] 
       
1.08
1.08
xa[3]=xa[2]+0.2*f(0.4,xa[2]);xa[3] 
       
1.2464
1.2464
xa[4]=xa[3]+0.2*f(0.6,xa[3]);xa[4] 
       
1.515968
1.515968
xa[5]=xa[4]+0.2*f(0.8,xa[4]);xa[5] 
       
1.91852288
1.91852288
xa 
       
(1.0, 1.0, 1.08, 1.2464, 1.515968, 1.91852288)
(1.0, 1.0, 1.08, 1.2464, 1.515968, 1.91852288)

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') cp5=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(xs2(t,0,1)),(t,0,1),thickness=2,color='red') (cp5+saprox0+saprox1+sol1).show(ymin = 0, ymax = 3,figsize=5) 
       

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.

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] 
       
1.04
1.04
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] 
       
1.165664
1.165664
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] 
       
1.3926255872
1.3926255872
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] 
       
1.75056237505
1.75056237505
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] 
       
2.28967260055
2.28967260055

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

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

Ejercicio 3. Aplica los métodos de Euler y Euler modificado con siete pasos para obtener el valor de la solución de la ecuación $x'(t)=-(t+1)(x^3(t)-3x(t)+1)$

a) en $t=1$ determinada por la condición inicial $x(0)=2$.

b) en $t=2$ determinada por la condición inicial $x(1)=2$

c) Dibujar el campo de pendientes

 
       
 
       

Problema. Llamaremos $a$ a la última letra de tu DNI y $b$ a la penúltima. Vamos a suponer que la población de atún rojo sigue el modelo de Malthus $p'(t)=(1+a)p(t)/10$, con $p(t)$ el número de atunes y el tiempo se mide en años. Sometemos a dicha población a una pesca estacional de $400(b+1)(1+cos(t/(2\pi)))$ unidades. Es decir, el modelo queda

$p'(t)=(1+a)p(t)/10-400(b+1)(1+\cos(t/(2\pi))).$

a) Si en el año de comienzo del estudio hay $12121$ atunes, ¿cuántos atunes se estiman a final de año?

b) ¿Se extinguirán los atunes?