%Ejemplo de clase (t^2y + y)dy=(ty^2 - t)dt
% - Introduzca la ecuación diferencial : 'Dy=(t*(y^2)-t)/((t^2)*y + y)'
% - Introduzca la condición y(a)=b : 'y(0.25)=1.25'
% - Introduzca la función de trabajo : (t*(y^2)-t)/((t^2)*y + y)
% - Introduzca la condición inicial : 1.25
% - Introduzca el valor de a : 0.25
% - Introduzca el valor de b : 0.30
% - Introduzca el tamaño de paso h : 0.01
%la soca con soluciones implicitas
%cuando la la ecuacion diferencial tiene solucion implicita quitar lineas
%33,34,103 y quitar linea 121 y poner 122
fprintf('\n');
clear all
fprintf(' -------------------------------------\n')
fprintf(' MÉTODO DE RUNGE-KUTTA DE ORDEN CUATRO\n')
fprintf(' -------------------------------------\n')
fprintf('\n');
syms t y
d=input(' - Introduzca la ecuación diferencial en terminos de (t,y) : ');
n=input(' - Introduzca la condición y(a)=b : ');
f1=input(' - Introduzca la función de trabajo : ');
ya=input(' - Introduzca la condición inicial : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tamaño de paso h : ');
fprintf('\n\n');
fprintf(' - La solución de la ecuación diferencial es : \n\n\n');
m = dsolve(d,n,'t');
pretty(m);
fprintf('\n\n\n');
f=f1;
w(1)=ya;
ti=a;
n=(b-a)/h;
q(1)=a;
v(1)=a;
d=0;
c=0;
g=0;
e=1;
fprintf('- w0 = %1.8f ',double(ya));
fprintf('\n\n');
for i=0:(n)
fprintf('---------------');
fprintf('\n');
fprintf('- Iteración: %1.0f\n',e);
fprintf('---------------');
fprintf('\n\n');
k1=h*subs(f,{t,y},{ti,w(i+1)});
fprintf('- K1 = h * f(t%1.0f,w%1.0f)',i,i);
fprintf('\n');
fprintf('- K1 = %1.15f * f(%1.15f,w%1.0f)',h,ti,i);
fprintf('\n');
fprintf('- K1 = %2.15f',double(k1))
fprintf('\n\n');
k2=h*subs(f,{t,y},{(ti+(h/2)),(w(i+1)+(k1/2))});
fprintf('- K2 = h * f(t%1.0f + h/2 , w%1.0f + K1/2)',i,i);
fprintf('\n');
fprintf('- K2 = %1.15f * f(%1.15f + %1.15f , w%1.0f + %1.15f)',double(h),double(ti),double(h/2),double(i),double((k1)/2));
fprintf('\n');
fprintf('- K2 = %2.8f',double(k2))
fprintf('\n\n');
k3=h*subs(f,{t,y},{(ti+(h/2)),(w(i+1)+(k2/2))});
fprintf('- K3 = h * f(t%1.0f + h/2 , w%1.0f + K2/2)',i,i);
fprintf('\n');
fprintf('- K3 = %1.9f * f(%1.15f + %1.15f , w%1.0f + %1.15f)',double(h),double(ti),double(h/2),double(i),double((k2)/2));
fprintf('\n');
fprintf('- K3 = %2.8f',double(k3))
fprintf('\n\n');
k4=h*subs(f,{t,y},{(ti+h),(w(i+1)+k3)});
fprintf('- K4 = h * f(t%1.0f + h , w%1.0f + K3)',i,i);
fprintf('\n');
fprintf('- K4 = %1.15f * f(%1.15f + %1.15f , w%1.0f + %1.15f)',double(h),double(ti),double(h),double(i),double(k3));
fprintf('\n');
fprintf('- K4 = %2.8f',double(k4))
fprintf('\n\n');
w(i+2)=w(i+1)+((1/6)*(k1+(2*k2)+(2*k3)+k4));
fprintf('- w%1.0f = w%1.0f + (1/6)*(K1+2K2+2K3+K4)',i+1,i)
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + (1/6)*(%1.15f+2*%1.15f+2*%1.15f+%1.15f)',double(i+1),double(i),double(k1),double(k2),double(k3),double(k4))
fprintf('\n');
fprintf('- w%1.0f = %2.8f',i+1,double(w(i+2)))
fprintf('\n\n');
ti=ti+h;
e=e+1;
end
fprintf('\n\n');
%Este for obtiene y guarda todos los valores de t
%También se utiliza para evaluar la ecuación diferencial
for p=a:h:b
d=1+d;
q(d)=p;
v(d)=subs(m,p);
end
%Este for se usa para contabilizar las iteraciones
for s=c:1:(d-1)
g=1+g;
k(g)=(g-1);
end
k3=k(end);
%Presentación de los datos
fprintf(' i ti wi+1 y(t)');
fprintf('\n\n');
for k1=0:k3
k2=k1+1;
fprintf('\n');
fprintf(' %1.0f %10.15f %10.15f %10.15f',k(k2),q(k2),w(k2),v(k2));
%fprintf(' %1.0f %10.15f %10.15f ',k(k2),q(k2),w(k2));
fprintf('\n');
end
fprintf('\n');
____________________________________________________________________________
%MÉTODO DE HEUN
%Ejemplo de clase y'=x*sen(2*x) + y*tan(x)
% - Introduzca la ecuación diferencial : 'Dy=x*sen(2*x) + y*tan(x)'
% - Introduzca la condición y(a)=b : 'y(2)=4'
% - Introduzca la función de trabajo : x*sen(2*x) + y*tan(x)
% - Introduzca la condición inicial : 4
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 0.2
% - Introduzca el tamaño de paso h : 0.05
%cuando la la ecuacion diferencial tiene solucion implicita quitar lineas
%35,36,55 y quitar linea 99 y poner 100
fprintf('\n');
clear all
fprintf(' --------------\n')
fprintf(' MÉTODO DE HEUN\n')
fprintf(' --------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuación diferencial : ');
n=input(' - Introduzca la condición y(a)=b : ');
f1=input(' - Introduzca la función de trabajo : ');
ya=input(' - Introduzca la condición inicial : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tamaño de paso h : ');
fprintf('\n\n');
fprintf(' - La solución de la ecuación diferencial es : \n\n\n');
m = dsolve(d,n,'x');
pretty(m);
fprintf('\n\n\n');
%Condiciones para el funcionamiento de los lazos FOR
f=f1;
w(1)=ya;
i=0;
t(1)=a;
v(1)=a;
d=0;
c=0;
g=0;
%Este for obtiene y guarda todos los valores de t
%También se utiliza para evaluar la ecuación diferencial
for p=a:h:b
d=1+d;
t(d)=p;
v(d)=subs(m,p);
end
%Este for se usa para contabilizar las iteraciones
for s=c:1:(d-1)
g=1+g;
k(g)=(g-1);
end
k3=k(end);
%Este for obtiene los valores aproximados de solución
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n');
fprintf(' FÓRMULAS DE CADA ITERACIÓN');
fprintf('\n');
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n\n');
fprintf('- w0 = %1.5f ',double(ya));
fprintf('\n');
for j=a:h:(b-h)
i=1+i;
w(i+1)=w(i)+((h/4)*(subs(f,{x,y},{t(i),w(i)})))+(((3/4)*h)*(subs(f,{x,y},{(t(i)+((2/3)*h)),(w(i)+(((2/3)*h)*(subs(f,{x,y},{t(i),w(i)}))))})));
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + h/4 f(t%1.0f,w%1.0f) + 3/4 h f(t%1.0f + 2/3 h,w%1.0f + 2/3 h f(t%1.0f,w%1.0f))',i,i-1,i-1,i-1,i-1,i-1,i-1,i-1);
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + %1.5f f(%1.15f,w%1.0f) + %1.5f f(%1.15f + %1.5f,w%1.0f + %1.5f f(%1.15f,w%1.0f))',i,i-1,h/4,t(i),i-1,(3/4)*h,t(i),(2/3)*h,i-1,(2/3)*h,t(i),i-1);
fprintf('\n');
end
fprintf('\n');
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n');
%Presentación de los datos
fprintf('\n\n');
fprintf(' i ti wi y(t)');
fprintf('\n\n');
for k1=0:k3
k2=k1+1;
fprintf('\n');
fprintf(' %1.0f %10.15f %10.15f %10.15f',k(k2),t(k2),w(k2),v(k2));
%fprintf(' %1.0f %10.15f %10.15f ',k(k2),t(k2),w(k2));
fprintf('\n');
end
fprintf('\n');
No hay comentarios:
Publicar un comentario