TAREA FINAL DE MÉTODOS NUMÉRICOS.
20%
Diseñar un programa en Matlab para resolver un problema de valor inicial de tercer orden
usando el método de Runge-Kutta de orden cuatro.
x ´ (t )=fx ( t , x , y , z ) , x ( ti )=xi
'
y (t )=fy ( t , x , y , z ) , y ( ti )= yi
'
z ( t )=fz ( t , x , y , z ) , z ( ti )=zi
En el intervalo: ti ≤ t ≤ tf con n particiones.
'' ' '' '
w ( t ) + p ( t ) w +q ( t ) w ( t ) +r ( t ) w ( t )=f (t)
' ''
w (ti )=a , w ( ti )=b , w ( ti )=c
x (t )=w ( t )
y ( t ) =w ' ( t )
z (t )=w ' ' ( t )
Ecuaciones:
x ' ( t )= y ( t )
y ' ( t )=z ( t )
'
z ( t )=f (t )−r ( t ) x ( t )−q ( t ) y ( t )− p ( t ) z (t )
Programa Runge-kutta 4 tercer orden
clear all
%Programa para resolver un problema de valor inicial de tercer orden
%por el método Runge-Kutta de orden cuatro
%La formulación del problema es:
%x'=fx(t,x,y,z);x(ti)=a;
%y'=fy(t,x,y,z);y(ti)=b;
%z'=fz(t,x,y,z);z(ti)=c;
%en e
l intervalo:[ti,tf] con n iteraciones.
% Entrada de la información
syms t x y z
n=input('Ingrese el número de iteraciones: ');
ti=input('Ingrese ti=');
tf=input('Ingrese tf=');
xi=input('Ingrese xi=');
yi=input('Ingrese yi=');
zi=input('Ingrese zi=');
f1=input('Ingrese la funcion fx(t,x,y,z)=');
f2=input('Ingrese la funcion fy(t,x,y,z)=');
f3=input('Ingrese la funcion fz(t,x,y,z)=');
h=(tf-ti)/n;
T=ti:h:tf;
X(1)=xi;
Y(1)=yi;
Z(1)=zi;
for k=1:n
t=T(k);
x=X(k);
y=Y(k);
z=Z(k);
M1x(k)=eval(f1);
M1y(k)=eval(f2);
M1z(k)=eval(f3);
t=T(k)+h/2;
x=X(k)+(h/2)*M1x(k);
y=Y(k)+(h/2)*M1y(k);
z=Z(k)+(h/2)*M1z(k);
M2x(k)=eval(f1);
M2y(k)=eval(f2);
M2z(k)=eval(f3);
x=X(k)+(h/2)*M2x(k);
y=Y(k)+(h/2)*M2y(k);
z=Z(k)+(h/2)*M2z(k);
M3x(k)=eval(f1);
M3y(k)=eval(f2);
M3z(k)=eval(f3);
t=T(k)+h;
x=X(k)+h*M3x(k);
y=Y(k)+h*M3y(k);
z=Z(k)+h*M3z(k);
M4x(k)=eval(f1);
M4y(k)=eval(f2);
M4z(k)=eval(f3);
Mx(k)=(1/6)*(M1x(k)+2*M2x(k)+2*M3x(k)+M4x(k));
My(k)=(1/6)*(M1y(k)+2*M2y(k)+2*M3y(k)+M4y(k));
Mz(k)=(1/6)*(M1z(k)+2*M2z(k)+2*M3z(k)+M4z(k));
X(k+1)=X(k)+h*Mx(k);
Y(k+1)=Y(k)+h*My(k);
Z(k+1)=Z(k)+h*Mz(k);
end
[T' X' Y' Z']
plot(T,X);
hold on
plot(T,Y,'c--');
plot(T,Z,'k--');
grid on
hold off
Ejercicio Equipo 3
'' ' '' ' −2 t
w ( t ) +4 w +8 w ( t )+ 8 w ( t )=t e
' ''
w ( 0 )=1 , w ( 0 )=1 , w ( 0 )=0
0 ≤ t ≤ 2, n=20
Entrada de información.
n= 20
ti= 0
tf= 2
xi= 1
yi= 1
zi= 0
fx= y
fy= z
fz= t*exp(-2*t)-4*z-8*y-8*x
Cálculos
h=(tf-ti)/n
T=ti:h:tf
X(k+1)=X(k)+h*Mx(k)
Y(k+1)=X(k)+h*My(k)
Z(k+1)=Z(k)+h*Mz(k)
Resultados.
[T’ X’ Y’ Z’]
0 1.0000 1.0000 0
0.1000 1.0976 0.9291 -1.3319
0.2000 1.1822 0.7493 -2.1921
0.3000 1.2453 0.5034 -2.6701
0.4000 1.2818 0.2253 -2.8461
0.5000 1.2902 -0.0582 -2.7910
0.6000 1.2707 -0.3272 -2.5669
0.7000 1.2257 -0.5677 -2.2272
0.8000 1.1584 -0.7704 -1.8170
0.9000 1.0730 -0.9300 -1.3734
1.0000 0.9739 -1.0449 -0.9261
1.1000 0.8655 -1.1159 -0.4984
1.2000 0.7521 -1.1458 -0.1071
1.3000 0.6376 -1.1390 0.2360
1.4000 0.5254 -1.1005 0.5239
1.5000 0.4184 -1.0361 0.7536
1.6000 0.3188 -0.9517 0.9249
1.7000 0.2285 -0.8530 1.0405
1.8000 0.1485 -0.7453 1.1048
1.9000 0.0796 -0.6336 1.1235
2.0000 0.0218 -0.5219 1.1032
Gráficas de las funciones en la misma figura.
Calculando la solución con dsolve
Podemos ver que las gráficas son equivalentes y que se superpone a la solución numérica
como se observa en la última gráfica.
