Sdgdf
Capítulo 2. Ecuaciones diferenciales ordinarias (EDOs). 2.1 Resolución de una ecuación diferencial ordinaria.
Vamos a resolver numéricamente la siguiente ecuación diferencial:
dy =y 1−y dt
con condición inicial y(0) = 0.2. Para empezar vamos a utilizar el método de Euler explícito. Este algoritmo, que estudiaréis en teoría con más detalle, se basa en la aproximación:
y tdt ≈y t dt ∙ y ' t =y tdt ∙ y t1−y t
para un valor de dt pequeño. Para empezar, tomamos un conjunto de instantes de tiempo t n=t 0 n dt y definimos y 0=y 0 . A continuación, definimos el valor de los puntos y n de forma recursiva, mediante el esquema numérico:
y n1=y n dt ∙ y n 1−y n Lo implementamos del siguiente modo tsteps=1000; % Subdivisiones del intervalo t_f=10; % valor final de tiempo dt=t_f/tsteps; % el paso de tiempo t=0:dt:t_f; % definimos los diferentes tiempos y=zeros(1,tsteps+1);% inicializamos y y(1)=0.2; % definimos dato inicial for n=1:tsteps y(n+1)=y(n)+dt*y(n)*(1y(n)); end %metodo de Euler explicito Dibujamos los resultados plot(t,y);Matlab tiene implementados varios métodos numéricos que podemos aplicar directamente. Crea un fichero de Matlab (File/New/Mfile) que contenga las siguientes líneas: function deriv=edo(t,y) deriv=y*(1y); Guarda el archivo con el nombre edo.m y establece el directorio donde lo has guardado como directorio de trabajo. Vuelve a la línea de comandos y teclea:
t_0=0; t_f=10; y0=0.2; [t,y]=ode23('edo',[t_0,t_f],y0); plot(t,y)La function Matlab ode23 es una implementación de un método RungeKutta de orden (2,3). En este caso no es necesario poner el paso de tiempo porque la misma función ode23 lo va reduciendo sobre la marcha si es necesario. Otra function Matlab es ode45 basado en un método explícito RungeKutta de orden (4,5). En general, ode45 es la mejor función para aplicar en un primer momento a muchos de los problemas:[t,y]=ode45('edo',[t_0,t_f],y0); plot(t,y) En este caso, es fácil comprobar que la función:
y t=
et 4et
es la solución exacta del problema de valor inicial anterior. Podemos usar esta referencia para comprobar qué tal precisión hemos obtenido con los distintos métodos. y_exacta=exp(t)./(4+exp(t)); plot(t,abs(yy_exacta)); %repetir cuando y se ha calculado con % distintos metodos
2.2 Errores.Vamos a estudiar la precisión que obtenemos en función del paso de tiempo utilizado. Para empezar creamos una función en matlab que devuelva sólo el valor aproximado de y(1) tomando como argumento el número de pasos temporales: function y=euler1(tsteps) t_f=1; % valor final de tiempo dt=t_f/tsteps; % el paso de tiempo t=0:dt:t_f; % definimos los diferentes tiempos y=0.2; % definimos dato inicialfor n=1:tsteps y=y+dt*y*(1y); end %metodo de Euler explicito Ahora podemos hacer una gráfica de la precisión obtenida (comparamos con el valor exacto,
e ≈0,404609675 ) en función del paso temporal dt. 4e
yexacto=exp(1)/(4+exp(1)); pasos=[];errores=[]; %vectores con los pasos temporales y los %errores for k=[4 10 30 100 300 1000] pasos=[pasos, 1/k]; %engordamos los vectores sobre la marchaerrores=[errores, abs(euler1(k)yexacto)]; end plot(pasos, errores)
En la siguiente gráfica a la izquierda puedes ver el resultado, mientras que a la derecha se muestra la gráfica correspondiente a un método de RungeKutta de orden 2 (el método del punto medio). Como puedes ver, la convergencia del método de Euler es lineal mientras que la del método del punto medio es cuadrática (la gráfica es aproximadamente una parábola).
En general, para estimar el orden de convergencia de un método, se escribe el error como una potencia del paso de tiempo: y se toma logaritmos para obtener una ecuación lineal:
error=C ∙ dt
log error= ∙ logdt logC
Al aplicar este procedimiento al error del método de Euler, se obtiene una recta (a la derecha, en ...
Regístrate para leer el documento completo.