Objetivo General:
* Aplicar los conocimientos básicos del cálculo, utilizando el lenguaje de programación Matlab.
Objetivos Específicos:
* Aplicar el algoritmo necesario para resolver ecuaciones diferenciales ordinarias (EDO), a través de una pequeña aplicación desarrollada en Matlab.
* Aplicar los algoritmos necesarios para resolver EDO, utilizando métodos numéricos, y en este caso particular, el método de Euler y Euler mejorado, a través de aplicaciones en Matlab.
Introducción
Las ecuaciones diferenciales aparecen naturalmente al modelar situaciones físicas en las ciencias naturales, ingeniería, y otras disciplinas, donde hay envueltas razones de cambio de una ó varias funciones desconocidas con respecto a una ó varias variables independientes. Estos modelos varían entre los más sencillos que envuelven una sola ecuación diferencial para una función desconocida, hasta otros más complejos que envuelven sistemas de ecuaciones diferenciales acopladas para varias funciones desconocidas. Por ejemplo, la ley de enfriamiento de Newton y las leyes mecánicas que rigen el movimiento de los cuerpos, al ponerse en términos matemáticos dan lugar a ecuaciones diferenciales. Usualmente estas ecuaciones están acompañadas de una condición adicional que especifica el estado del sistema en un tiempo o posición inicial. Esto se conoce como la condición inicial y junto con la ecuación diferencial forman lo que se conoce como el problema de valor inicial. Por lo general, la solución exacta de un problema de valor inicial es imposible ó difícil de obtener en forma analítica. Por tal razón los métodos numéricos se utilizan para aproximar dichas soluciones. En este caso utilizaremos los métodos de Euler y Euler mejorado.
Ecuaciones Diferenciales Numéricas
METODO DE EULER
Considere el problema de valor inicial para la función (desconocida) y(t) descrito por:
Defina para n³0 los siguientes cantidades:
Para cualquier j³0, tenemos por el Teorema de Taylor que podemos escribir:
Codigo:
%metodo de euler para 1 ec dif
function [t,y]=euler(f,I,y0,h)
n=(I(2)-I(1))/h;
y(1)=y0;
t=I(1):h:I(2);
for i=1:n
y(i+1)=y(i)+h*feval(f,t(i),y(i));
end
t=t';
y=y';
plot(t,y,'r')
title('METODO DE EULER PARA UNA EDO')
xlabel('t')
ylabel('y')
%por ejemplo:Aplicar el metodo de Euler para aproximar la ecuacion:
% y'(t)=y + 3.t^2 + 2, y(0) = 1 en el intervalo [0,2] con h = 0,2.
%ejecucion:
%[t,y]=euler('epa',[0,2],1,0.2)
%poner en otro editor de texto el sgte programa con nombre epa.m donde defines
%la edo a resolver
function f=epa(t,u)
f(1)=u(1)+3*t^2+2;
MÉTODO DE HEUN (GRAL)
Ecuación:
Para propositos de hacer cálculos es mejor escribir esta fórmula como:
Código:
clc
% METODO DE HEUN PARA SISTEMA DE ECUACIONES DIFERENCIALES.
function [t,y]=heunbc(f,I,y0,n)
h=intp(I,n);% LLAMANDO A intp.m
t=partir(I,n);% LLAMANDO A partir.m
y(1,:)=y0;
h2=h/2;
for i=1:n
k1=feval(f,t(i), y(i,:));
y(i+1,:)=y(i,:)+h2*(k1+feval(f,t(i)+h,y(i,:)+h*k1));
end
t=t';
%ejemplo:
%Resolver el sistema de ecuaciones diferenciales dado por
%u'(1)=u2 0<=t<=2 u1(0) = 1
%u'(2)=-u1-2e^t + 1 0<=t<=2 u2(0) = 0
%u'(3)=-u1-e^t + 1 0<=t<=2 u3(0) = 1
%%EJECUCION: [t,y]=heunbc('veo',[0,2],[1 0 1],0.2)
%esto copia en un editor y guarda con el mismmo nombre
function f=veo(t,u)
f(1)=u(2);
f(2)=-u(1)-2*exp(t)+1;
f(3)=-u(1)-exp(t)+1;
% EN "VEO" SE DEFINe LA EDO Q DESEAS RESOLVER
%esto copia en otro editor y guarda con el mismo nombrre
function y=intp(I,noh)
y=(I(2)-I(1))/(noh);
%esto copia en otro editor y guarda con el mismo nombrre
function y=partir(I,n)
h=intp(I,n);
for i=1:n+1
X(i)=I(1)+(i-1)*h;
end
y=X;
%Ahora si puedes ejecutar el programa en el comando window
%pon la ejecucion q esta en la parte superior.comenta... Ejemplo:
Consideremos aqui las ecuaciones diferenciales que se obtienen de las leyes de Newton aplicadas al problema de dos cuerpos. Suponemos aqui que uno de los cuerpos es mucho más masivo que el otro de modo que su movimiento es descartable, e.g., la tierra y un satélite. Suponemos también que el movimiento es en un plano. Como la fuerza gravitacional es inversamente proporcional a la distancia entre los cuerpos, tenemos tomando todas las constantes envueltas como uno, que la posición (x(t),y(t)) del cuerpo pequeño esta dada por el sistema de ecuaciones diferenciales:
Tomamos como condiciones iniciales . Debido a que este es un sistema de orden dos, tenemos que hacer la sustitución:
lo cual tranforma el sistema de arriba al siguiente sistema de orden uno:
clear all
disp('METODO DE EULER MODIFICADO')
clc
syms x
syms y
f=inline(input('ingrese la derivada:','s'));
x=input('ingrese el valor de x:');
y=input('ingrese el valor de y:');
h=input('ingrese el valor de h:');
n=input('ingrese numero de iteraciones:');
clc
disp('x(n) y´(n) hy´(n) y(n+1),p hy´(n+1),p y(n+1),c');
for i=1:n
s=h+x;
y1=feval(f,x,y);
hy1=h*y1;
y2=y+hy1;
y3=feval(f,s,y2);
hy2=y3*h;
yn=y+((hy1+hy2)/2);
fprintf('n%0.1f %0.4f %0.4f %0.4f %0.4f %0.4f',x,y,hy1,y2,hy2,yn);
y=yn;
x=x+h;
end
METODO DE LA SECANTE
clc
syms x
f=inline(input('ingrese la funcion :','s'));
x0=input('ingrese intervalo inferior:');
x1=input('ingrese intervalo superior:');
tol=input('ingrese la tolerancia:');
fx0=feval(f,x0);
fx1=feval(f,x1);
n=1;
if(abs(fx0))<(abs(fx1)) fx2=900; disp('n x0 x1 x2 f(x0) f(x1) f(x2) ') while (abs(fx2))>tol
x2=x1-((x1-x0)/(fx1-fx0))*fx1;
fx2=feval(f,x2);
fprintf('n%d %0.4f %0.4f %0.4f %0.4f %0.4f %0.4f',n,x0,x1,x2,fx0,fx1,fx2)
x0=x1;
x1=x2;
fx0=feval(f,x0);
fx1=feval(f,x1);
n=n+1;
end
else
disp('NO SE PUEDE REALIZAR POR ESTE METODO');
disp('|F(X0)| TIENE QUE SER MENOR A |F(X1)|');
end
METODO DE RUUGE-KUTTA
Uno de los métodos más utilizados para resolver numéricamente problemas de ecuaciones diferenciales ordinarias con condiciones iniciales es el método de Runge-Kutta de cuarto orden, el cual proporciona un pequeño margen de error con respecto a la solución real del problema y es fácilmente programable en un software para realizar las iteraciones necesarias.
clc
disp('METODO DE RUUGE-KUTTA')
syms x
syms y
f=inline(input('ingrese la derivada:','s'));
x=input('ingrese el valor de x:');
y=input('ingrese el valor de y:');
h=input('ingrese el valor de h:');
n=input('ingrese numero de iteraciones:');
clc
disp('x(n) k1 k2 k3 k4 ');
for i=1:n
k1=h*(feval(f,x,y));
z=x+(1/2*(h));
w=y+(1/2*(k1));
k2=h*(feval(f,z,w));
w=y+(1/2*(k2));
k3=h*(feval(f,z,w));
z=x+h;
w=y+k3;
k4=h*(feval(f,z,w));
fprintf('n%0.1f %0.4f %0.4f %0.4f %0.4f %0.4f',x,y,k1,k2,k3,k4);
x=x+h;
y=y+((1/6)*(k1+2*k2+2*k3+k4));
end
f=inline(input('ingrese la funcion:','s'));
a=input('ingrese intervalo inferior:');
b=input('ingrese intervalo superior:');
tol=input('ingrese la tolerancia:');
fa=feval(f,a);
fb=feval(f,b);
if (sign(fa)==sign(fb))
disp('no se garantiza el teorma de bolzano:');
break;
end
fc=10;
n=1;
disp('n an bn cn f(c)n');
while(n~=8)
c=(a+b)/2;
fc=feval(f,c);
fprintf('n%d %0.4f %0.4f %0.4f %0.4f',n,a,b,c,fc);
if (fc==0)
fprintf('%f es un cero de f',c);
break;
end
if (sign(fc)==sign(fa))
s=c;
t=b;
else
if(sign(fc)~=sign(fb))
s=a;
t=c;
end
end
n=n+1;
a=s;
b=t;
fa=feval(f,a);
fb=feval(f,b);
end
POSDATA: no se les olvide agradecer...............
. Debido a que este es un sistema de orden dos, tenemos que hacer la sustitución: 
