GT_MM_1314_Ecuaciones_Diferenciales

1395 days ago by trinidad

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

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+plot(sol0(t,0.3),(t,-2,2),thickness=4,color='green')+plot(sol0(t,0),(t,-2,2),thickness=4,color='blue')+plot(sol0(t,-0.5),(t,-2,2),thickness=4,color='red')).show(ymin = -2, ymax = 2) 
       

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 # Incremento de población por unidad de tiempo sf1=plot_slope_field(Tc*(1-x/K)*x, (t,0,5), (x,0,K*2));sf1 
       

Vamos a dibujar dos 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)) sf1+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') 
       

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.


 
       
# Solución plot_slope_field(cos(t)*x, (t,0,2*pi), (x,-2,2))+plot(exp(sin(t))/2,(t,0,2*pi),color='red')+plot(-1-sin(t),(t,0,2*pi)) 
       

2. Ecuaciones autónomas

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

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') cp1=plot_slope_field(f(x), (t,-2,2), (x,-1,2));cp1 
       

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 
       

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

var('c') solve(h==t+c,x) 
       

Guardamos la solución en una variable.

xs(t,c)= e^(c + t)/(e^(c + t) - 1);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) 
       

Dibujamos la función en el campo de pendientes:

sol03=plot(xs(t,log(3/2)),(t,-2,2),thickness=4,color='green') (cp1+sol03).show(ymin = -1, ymax = 2) 
       

Como la función que hemos definido como solución tiene una asíntota, observamos dos ramas. De nuevo podemos utilizar solve para obtener la posición de la asíntota y dibujar una única rama:

xs(t,log(3/2)) 
       
solve(e^(t + log(3/2)) - 1==0,t) 
       
sol03=plot(xs(t,log(3/2)),(t,log(2/3),2),thickness=4,color='green') (cp1+sol03).show(ymin = -1, ymax = 2) 
       

Vamos a intertar resolver ahora el problema de valor inicial $x'(t)=x(t)(1-x(t))$, $x(0)=1/2$. Seguimos el mismo procedimiento:

solve(xs(0,c)==1/2,c) 
       

Sage nos devuelve un número complejo. Recordemos que en clase lo que hacíamos era que al despejar teníamos dos posibilidades y descartábamos una al imponer la condición inicial. Sage no necesita hacer eso porque trabaja con números complejos. Podemos sustituir el valor de c que hemos obtenido en la solución y comprobar que la solución que tenemos es precisamente la que queremos:

simplify(xs(t,I*pi)) 
       
sol012=plot(simplify(xs(t,I*pi)),(t,-2,2),thickness=4,color='red') (cp1+sol03+sol012).show(ymin = -1, ymax = 2) 
       

Ejercicio 2. 

a) Calcula los puntos críticos de $x'=x(1-x)$.

b) Calcula las soluciones de la ecuación $x'=x(1-x)$ determinadas por las condiciones iniciales $x(0)=-0.5$ y $x(2)=0.5$.

c) Añade a la gráfica anterior (del campo de pendientes con las soluciones) las soluciones obtenidas en los apartados a) y b).

 
       

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') cp2=plot_slope_field(f(x)*g(t), (t,0,2), (x,-pi,pi));cp1 
       

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

ix=integrate(1/f(x),x).simplify_full() it=integrate(1/t,t).simplify_full() ix,it 
       

De nuevo utilizamos solve para despejar la $x$.

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

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

solve(xs(1,c)==3,c) 
       

Mostramos la solución:

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

Calculamos hasta donde está definida y la pintamos.

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

4. Ecuaciones lineales

Resolviendo $x'(t)=tx(t)+t$, tenemos

\[\left(e^{-A(t)}x(t) \right)'=te^{-A(t)}\] e integrando y despejando $x(t)$, \[x(t)=x(0)e^{A(t)}+e^{A(t)}\int_{t_0}^t x e^{-A(s)}\,ds. \]

s=var('s') A(t,t0)=integrate(s,s,t0,t) xs(t,t0,x0)=x0*exp(A(t,t0))+exp(A(t,t0))*integrate(s*exp(-A(s,t0)),s,t0,t);show(simplify(xs)) 
       
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=4,color='red') solb=plot(simplify(xs(t,0,0)),(t,-2,2),thickness=4,color='green') (cp3+sola+solb).show(ymin = -2, ymax = 2) 
       

Problema 1. 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? 

 
       

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

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 
       

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=4,color='green') saprox1=line(zip(vector(RDF,range(0,6))*0.2,xa),thickness=4,color='blue') sol1=plot(simplify(xs(t,0,1)),(t,0,1),thickness=4,color='red') (cp4+saprox0+saprox1+sol1).show(ymin = 0, ymax = 3) 
       

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] 
       
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=4,color='purple') (cp4+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 
       

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') (cp4+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 anteriores con tres pasos para obtener el valor de la solución de la ecuación $x'(t)=(t^2+1)(x^3(t)+2x(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

 
       

Rellenar el cuadro de texto, imprimir a pdf y enviar la tarea al campus virtual.

Alumnos: 

 
       

Curvas geodésicas

A continuación vamos a mostrar una pequeña aplicación en Sage que dibuja las curvas geodésicas en una superficie. Está tomado de http://wiki.sagemath.org/interact/geometry.

Representación de la superficie.

u, v, t = var('u v t') @interact def _(x = input_box(3*sin(u)*cos(v), 'x'), y = input_box(sin(u)*sin(v), 'y'), z = input_box(2*cos(u), 'z'), _int_u = input_grid(1, 2, default = [[0,pi]], label = 'u -interval'), _int_v = input_grid(1, 2, default = [[-pi,pi]], label = 'v -interval')): global F, Fu, Fv, func, S_plot, int_u, int_v int_u = _int_u[0] int_v = _int_v[0] F = vector([x, y, z]) S_plot = parametric_plot3d( F, (u, int_u[0], int_u[1]), (v, int_v[0], int_v[1])) S_plot.show(aspect_ratio = [1, 1, 1]) dFu = F.diff(u) dFv = F.diff(v) Fu = fast_float(dFu, u, v) Fv = fast_float(dFv, u, v) ufunc = function('ufunc', t) vfunc = function('vfunc', t) dFtt = F(u=ufunc, v=vfunc).diff(t, t) ec1 = dFtt.dot_product(dFu(u=ufunc, v=vfunc)) ec2 = dFtt.dot_product(dFv(u=ufunc, v=vfunc)) dv, ddv, du, ddu = var('dv, ddv, du, ddu') diffec1 = ec1.subs_expr(diff(ufunc, t) == du, diff(ufunc, t, t) == ddu, diff(vfunc, t) == dv, diff(vfunc, t, t) == ddv, ufunc == u, vfunc == v) diffec2 = ec2.subs_expr(diff(ufunc, t) == du, diff(ufunc, t, t) == ddu, diff(vfunc, t) == dv, diff(vfunc, t, t) == ddv, ufunc == u, vfunc == v) sols = solve([diffec1 == 0 , diffec2 == 0], ddu, ddv) ddu_rhs = (sols[0][0]).rhs().full_simplify() ddv_rhs = (sols[0][1]).rhs().full_simplify() ddu_ff = fast_float(ddu_rhs, du, dv, u, v) ddv_ff = fast_float(ddv_rhs, du, dv, u, v) def func(y,t): v = list(y) return [ddu_ff(*v), ddv_ff(*v), v[0], v[1]] 
       

Click to the left again to hide and once more to show the dynamic interactive window

Representación de la geodésica que pasa por un punto, con cierta dirección.

# second interact from scipy.integrate import odeint def fading_line3d(points, rgbcolor1, rgbcolor2, *args, **kwds): L = len(points) vcolor1 = vector(RDF, rgbcolor1) vcolor2 = vector(RDF, rgbcolor2) return sum(line3d(points[j:j+2], rgbcolor = tuple( ((L-j)/L)*vcolor1 + (j/L)*vcolor2 ), *args, **kwds) for j in srange(L-1)) steps = 100 @interact def _(u_0 = slider(int_u[0], int_u[1], (int_u[1] - int_u[0])/100, default = (int_u[0] + int_u[1])/2, label = 'u_0'), v_0 = slider(int_v[0], int_v[1], (int_v[1] - int_v[0])/100, default = (int_v[0] + int_v[1])/2, label = 'v_0'), V_u = slider(-10, 10, 1/10, default = 1, label = 'V_u'), V_v = slider(-10, 10, 1/10, default = 0, label = 'V_v'), int_s = slider(0, 10, 1/10, default = (int_u[1] - int_u[0])/2, label = 'geodesic interval'), sliding_color = checkbox(True,'change color along the geodesic')): du, dv, u, v = var('du dv u v') Point = [u_0, v_0] velocity = [V_u, V_v] Point = map(float, Point) velocity = map(float, velocity) geo2D_aux = odeint(func, y0 = [velocity[0], velocity[1], Point[0], Point[1]], t = srange(0, int_s, 0.01)) geo3D = [F(u=l,v=r) for [j, k, l, r] in geo2D_aux] if sliding_color: g_plot = fading_line3d(geo3D, rgbcolor1 = (1, 0, 0), rgbcolor2 = (0, 1, 0), thickness=4) else: g_plot = line3d(geo3D, rgbcolor=(0, 1, 0), thickness=4) P = F(u=Point[0], v=Point[1]) P_plot = point3d((P[0], P[1], P[2]), rgbcolor = (0, 0, 0), pointsize = 30) V = velocity[0] * Fu(u = Point[0], v = Point[1]) + \ velocity[1] * Fv(u= Point[0], v = Point[1]) V_plot = arrow3d(P, P + V, color = 'black') show(g_plot + S_plot + V_plot + P_plot,aspect_ratio = [1, 1, 1]) 
       

Click to the left again to hide and once more to show the dynamic interactive window