jueves, 22 de mayo de 2014
uso de los comandos ode, ejm ode 45
[t y]=ode45(@dy,[a:h:b],y(a))
en el editor se escribe
function dy=g(t,y)
dy=2*sin(2*t);
[t y]=ode45(@g,[1:1:5],2)
para un sistema de ecuaciones despejar x' y y' o las variables que tenga y luego
[t, xy]=ode23(@fsis,[a:h:b],[x(a);y(a)])
para el archivo .m
function p=fsis(t,w)
p=[coeficientes de las variables de x' dejando espacios ; igual con y']*w+[numeros y termnos de t]
metodo modificado de euler
f=input('ingrese la funcion de trabajo entre comillas');
g=input('ingrese la ED entre comillas');
e=input('ingrese la condicion inicial y(a)= entre comillas');
a=input('ingrese el valor a');
b=input('ingrese el valor b');
ya=input('ingrese la condicion inicial (el valor de y(a)) ');
h=input('ingrese h );
m=(b-a)/h;
T=a:h:b;
p=dsolve(g,e);
w(1)=ya;
for j==m+1
w(j+1)=w(j) + 0.5*h*subs(f,{t,y},{T(j),w(j)}) + 0.5*h*subs(f,{t,y},{T(j+1),w(j) + h*subs(f,{t,y},{T(j),w(j)})})
end
Y=double(subs(p,t,T));
fprintf(' T wi+1' Y(t));
fprintf('\n ');
E=[T' w' Y'];
si probar con programas mas abajo
%METODO DE NEVILLE
fprintf('\n');
clear all
clc
fprintf(' -----------------\n')
fprintf(' MÉTODO DE NEVILLE\n')
fprintf(' -----------------\n')
fprintf('\n');
syms x
res=input('- La Funcion le fue dada(Si=1,No=0)? : ');
if res==1
fun=input('- Introduzca la Funcion F(x) : ','s');
end
Xi=input('- Introduzca la cantidad para aproximar : ');
n=input('- Introduzca la cantidad de puntos dados : ');
fprintf('\n\n');
for i=0:(n-1),
fprintf('- Introduzca X%1.0f ',i);
X(i+1)=input (' = ');
if res==0
fprintf('- Introduzca F(X%1.0f) ',i);
FX(i+1)=input(' = ');
else
FX(i+1)=funcion(X(i+1),fun);
end
end
for i=1:n,
Q(i,1)=FX(i);
end
for i=2:n,
for j=i:n,
Q(j,i)=(((Xi-X(j-i+1))*Q(j,i-1))-((Xi-X(j))*Q(j-1,i-1)))/(X(j)-X(j-i+1));
end
end
fprintf('\n\n');
for i=2:n,
for j=i:n,
fprintf('- Q (%1.0f,%1.0f) = %3.8f\n ',j-1,i-1,Q(j,i));
fprintf('\n');
end
end
-----------------------------------------------------------------------------------------
%MÉTODO DEL PUNTO MEDIO
fprintf('\n');
clear all
clc
fprintf(' ----------------------\n')
fprintf(' MÉTODO DEL PUNTO MEDIO\n')
fprintf(' ----------------------\n')
fprintf('\n');
syms x y
disp('Wi+1=Wi+hF(Ti+h/2,Wi+h/2F(Ti,Wi)) la funcion de trabajo Dy/Dx=F(Ti,Wi)')
fprintf('\n');
d=input(' - Introduzca la ecuación diferencial encerradas en tildes y q Dy no este x y : ');
n=input(' - Introduzca la condición y(a)=b encerradas en tildes : ');
l=input(' - Introduzca la variable encerradas en tildes : ');
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,l);
pretty(m);
fprintf('\n\n\n\n');
fprintf(' i ti wi y(t)');
%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
%Este for obtiene los valores aproximados de solución
for j=a:h:(b-h)
i=1+i;
w(i+1)=w(i)+(h*(subs(f,{x,y},{(t(i)+h/2),w(i)+((h/2)*(((subs(f,{x,y},{t(i),w(i)})))))})));
end
%Presentación de los datos en forma de tabla utilizando matrices
[k' t' w' v']
----------------------------------------------------------------------------------------------------------------
%MÉTODO DE HEUN
fprintf('\n');
clear all
clc
fprintf(' --------------\n')
fprintf(' MÉTODO DE HEUN\n')
fprintf(' --------------\n')
fprintf('\n');
syms x y
disp('funcion de trabajo F(Ti,Wi)=Dy/Dx')
fprintf('\n');
disp('formula del metodo: Wi+1=Wi+h/4*F(Ti,Wi)+3h/4*F(Ti+2h/3, Wi+2h/3*F(Ti,Wi) )')
d=input(' - Introduzca la ecuación diferencial encerrados en tildes y q DY no est multiplpcd * y : ');
n=input(' - Introduzca la condición y(a)=b : ');
l=input(' - Introduzca la variable : ');
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,l);
pretty(m);
fprintf('\n\n\n\n');
fprintf(' i ti wi y(t)');
%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
%Este for obtiene los valores aproximados de solució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)}))))})));
end
%Presentación de los datos en forma de tabla utilizando matrices
[k' t' w' v']
_-----------------------------------------------------------------------------------------------------------
%MÉTODO DE RUNGE-KUTTA-FEHLBERG ORDEN 5
% - Introduzca la ecuación diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b : 'y(0)=0.5'
% - Introduzca la función de trabajo : y-(x^2)+1
% - Introduzca la condición inicial : 0.5
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tamaño de paso h : 0.1
clear all
disp('---------------------------------------------------------------------------');
disp('Método de RUNGE-KUTTA-FEHLBERG de orden 5');
disp('---------------------------------------------------------------------------');
fprintf('\n')
syms x y
d=input(' - Ingrese la ecuación diferencial en (x,y): ');
n=input(' - Ingrese la condición y(a)=b: ');
f1=input(' - Ingrese la función de trabajo: ');
ya=input(' - Ingrese el valor de la condicion inicial: ');
a=input(' - Ingrese el valor de a: ');
b=input(' - Ingrese el valor de b: ');
h=input(' - Ingrese el valor de h: ');
fprintf('\n');
fprintf(' - La solución de la ecuación diferencial es : \n\n\n');
m = dsolve(d,n,'x');
pretty(m);
fprintf('\n');
f=f1;
w(1)=ya;
i=0;
t(1)=a;
q(1)=a;
v(1)=a;
d=0;
c=0;
g=0;
e=1;
disp('---------------------------------------------------------------------------');
disp('Formulas para cada iteracion');
disp('---------------------------------------------------------------------------');
fprintf('\n')
fprintf('- w0 = %1.9f ',ya);
fprintf('\n');
for j=a:h:(b-h)
i=1+i;
t(i)=j;
fprintf('- Iteración: %1.0f\n',e);
fprintf('\n');
k1=h*subs(f,{x,y},{t(i),w(i)});
fprintf('- K1 = h * f(t%1.0f,w%1.0f)',i-1,i-1);
fprintf('\n');
fprintf('- K1 = %1.9f * f(%1.9f,w%1.0f)',h,t(i),i-1);
fprintf('\n');
fprintf('- K1 = %2.15f',double(k1))
fprintf('\n\n');
k2=h*subs(f,{x,y},{t(i)+h/4,w(i)+k1/4});
fprintf('- K2 = h * f(t%1.0f + h/4 , w%1.0f + K1/4)',i-1,i-1);
fprintf('\n');
fprintf('- K2 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f)',h,t(i),h/4,i-1,double(k1)/4);
fprintf('\n');
fprintf('- K2 = %2.15f',double(k2))
fprintf('\n\n');
k3=h*subs(f,{x,y},{t(i)+(3.*h/8),w(i)+(3.*k1/32)+(9.*k2/32)});
fprintf('- K3 = h * f(t%1.0f + (3/8)h , w%1.0f + 3K1/32 + 9K2/32)',i-1,i-1);
fprintf('\n');
fprintf('- K3 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f + %1.9f)',h,t(i),((3/8)*h),i-1,(3*(double(k1))/32),(9*(double(k2))/32));
fprintf('\n');
fprintf('- K3 = %2.15f',double(k3))
fprintf('\n\n');
k4=h*subs(f,{x,y},{t(i)+(12.*h/13),w(i)+(1932.*k1/2197)-(7200.*k2/2197)+(7296.*k3/2197)});
fprintf('- K4 = h * f(t%1.0f + (12/13)h , w%1.0f + 1932K1/2197 - 7200K2/2197 + 7296K3/2197)',i-1,i-1);
fprintf('\n');
fprintf('- K4 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f - %1.9f + %1.9f)',h,t(i),((12/13)*h),i-1,(1932*(double(k1))/2197),(7200*(double(k2))/2197),(7296*(double(k3))/2197));
fprintf('\n');
fprintf('- K4 = %2.15f',double(k4))
fprintf('\n\n');
k5=h*subs(f,{x,y},{t(i)+h,w(i)+(439.*k1/216)-(8.*k2)+(3680.*k3/513)-(845.*k4/4104)});
fprintf('- K5 = h * f(t%1.0f + h , w%1.0f + 439K1/216 - 8K2 + 3680K3/513 - 845K4/4104)',i-1,i-1);
fprintf('\n');
fprintf('- K5 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f - %1.9f + %1.9f - %1.9f)',h,t(i),h,i-1,(439*(double(k1))/216),(8*(double(k2))),(3680*(double(k3))/513),(845*(double(k3))/4104));
fprintf('\n');
fprintf('- K5 = %2.15f',double(k5))
fprintf('\n\n');
k6=h*subs(f,{x,y},{t(i)+h/2,w(i)-(8.*k1/27)+(2*k2)-(3544.*k3/2565)+(1859.*k4/4104)-(11.*k5/40)});
fprintf('- K6 = h * f(t%1.0f + h/2 , w%1.0f + 8K1/27 + 2K2 - 3544K3/2565 + 1859K4/4104 - 11K5/40)',i-1,i-1);
fprintf('\n');
fprintf('- K6 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f + %1.9f - %1.9f + %1.9f - %1.9f)',h,t(i),h,i-1,(8*(double(k1))/27),(2*(double(k2))),(3544*(double(k3))/2565),(1859*(double(k3))/4104),(11*(double(k3))/40));
fprintf('\n');
fprintf('- K6 = %2.15f',double(k6))
fprintf('\n\n');
w(1+i)=w(i)+(16.*k1/135)+(6656.*k3/12825)+(28561.*k4/56430)-(9.*k5/50)+(2.*k6/55);
fprintf('- w%1.0f = w%1.0f + (16/135) K1 + (6656/12825) K3 + (28561/56430) K4 - (9/50) K5 + (2/55) K6)',i,i-1)
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + %1.9f + %1.9f + %1.9f - %1.9f + %1.9f',i,i-1,(16/135)*double(k1),(6656/12825)*double(k3),(28561/56430)*double(k4),(9/50)*double(k5),(2/55)*double(k6))
fprintf('\n');
fprintf('- w%1.0f = %2.15f',i,double(w(i+1)))
fprintf('\n\n');
e=e+1;
end
fprintf('\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) Error');
fprintf('\n\n');
for k1=0:k3
k2=k1+1;
fprintf('\n');
fprintf(' %1.0f %10.9f %10.15f %10.15f %10.15f',k(k2),q(k2),w(k2),v(k2),abs(w(k2)-v(k2)));
fprintf('\n');
end
fprintf('\n');
----------------------------------------------------------------------------------------------------
_--------------------------------------------------------------------------------------------------
%MÉTODO DE RUNGE-KUTTA-FEHLBERG 4
% - Introduzca la ecuación diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b : 'y(0)=0.5'
% - Introduzca la función de trabajo : y-(x^2)+1
% - Introduzca la condición inicial : 0.5
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tamaño de paso h : 0.1
fprintf('\n');
clc
fprintf(' -----------------------------------------\n')
fprintf(' METODO DE RUNGE-KUTTA-FEHLBERG DE ORDEN 4\n')
fprintf(' -----------------------------------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuacion diferencial en (x,y) : ');
n=input(' - Introduzca la condicion y(a)=b : ');
f1=input(' - Introduzca la funcion de trabajo : ');
ya=input(' - Introduzca la condicion inicial w0 : ');
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 solucion de la ecuacion diferencial es : \n\n\n');
m = dsolve(d,n,'x');
pretty(m);
fprintf('\n\n\n');
f=f1;
w(1)=ya;
i=0;
t(1)=a;
q(1)=a;
v(1)=a;
d=0;
c=0;
g=0;
e=1;
fprintf('- w0 = %1.9f ',ya);
fprintf('\n\n');
for j=a:h:(b-h)
i=1+i;
t(i)=j;
fprintf(' ---------------');
fprintf('\n');
fprintf('- Iteracion: %1.0f\n',e);
fprintf(' ---------------');
fprintf('\n\n');
k1=h*subs(f,{x,y},{t(i),w(i)});
fprintf('- K1 = h * f(t%1.0f,w%1.0f)',i-1,i-1);
fprintf('\n');
fprintf('- K1 = %1.9e * f(%1.9e,w%1.0e)',h,t(i),i-1);
fprintf('\n');
fprintf('- K1 = %2.15e',(k1))
fprintf('\n\n');
k2=h*subs(f,{x,y},{t(i)+h/4,w(i)+k1/4});
fprintf('- K2 = h * f(t%1.0f + h/4 , w%1.0f + K1/4)',i-1,i-1);
fprintf('\n');
fprintf('- K2 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f)',h,t(i),h/4,i-1,(k1)/4);
fprintf('\n');
fprintf('- K2 = %2.15f',(k2))
fprintf('\n\n');
k3=h*subs(f,{x,y},{t(i)+(3.*h/8),w(i)+(3.*k1/32)+(9.*k2/32)});
fprintf('- K3 = h * f(t%1.0f + (3/8)h , w%1.0f + 3K1/32 + 9K2/32)',i-1,i-1);
fprintf('\n');
fprintf('- K3 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f + %1.9f)',h,t(i),((3/8)*h),i-1,(3*((k1))/32),(9*((k2))/32));
fprintf('\n');
fprintf('- K3 = %2.15f',(k3))
fprintf('\n\n');
k4=h*subs(f,{x,y},{t(i)+(12.*h/13),w(i)+(1932.*k1/2197)-(7200.*k2/2197)+(7296.*k3/2197)});
fprintf('- K4 = h * f(t%1.0f + (12/13)h , w%1.0f + 1932K1/2197 - 7200K2/2197 + 7296K3/2197)',i-1,i-1);
fprintf('\n');
fprintf('- K4 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f - %1.9f + %1.9f)',h,t(i),((12/13)*h),i-1,(1932*((k1))/2197),(7200*((k2))/2197),(7296*((k3))/2197));
fprintf('\n');
fprintf('- K4 = %2.15f',(k4))
fprintf('\n\n');
k5=h*subs(f,{x,y},{t(i)+h,w(i)+(439.*k1/216)-(8.*k2)+(3680.*k3/513)-(845.*k4/4104)});
fprintf('- K5 = h * f(t%1.0f + h , w%1.0f + 439K1/216 - 8K2 + 3680K3/513 - 845K4/4104)',i-1,i-1);
fprintf('\n');
fprintf('- K5 = %1.9f * f(%1.9f + %1.9f , w%1.0f + %1.9f - %1.9f + %1.9f - %1.9f)',h,t(i),h,i-1,(439*((k1))/216),(8*((k2))),(3680*((k3))/513),(845*((k3))/4104));
fprintf('\n');
fprintf('- K5 = %2.15f',(k5))
fprintf('\n\n');
w(1+i)=w(i)+(25*k1/216)+(1408*k3/2565)+(2197*k4/4104)-(k5/5);
fprintf('- w%1.0f = w%1.0f + (25/216) K1 + (1408/2565) K3 + (2197/4104) K4 - (1/5) K5)',i,i-1)
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + %1.9f + %1.9f + %1.9f - %1.9f',i,i-1,(25/216)*(k1),(1408/2565)*(k3),(2197/4104)*(k4),(1/5)*(k5))
fprintf('\n');
fprintf('- w%1.0f = %2.15f',i,(w(i+1)))
fprintf('\n\n');
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.15e %10.15e %10.15e',k(k2),q(k2),w(k2),v(k2));
fprintf('\n');
end
fprintf('\n');
--------------------------------------------------------------------------------------------
___________________________________________________
fprintf('\n');
clear all
fprintf(' --------------------------\n')
fprintf(' MÉTODO MODIFICADO DE EULER\n')
fprintf(' --------------------------\n')
fprintf('\n');
syms x y
disp('funcion de trabajo F(Ti,Wi)=Dy/DX')
fprintf('\n');
disp('formula del metodo: Wi+1=Wi+h/2*F(Ti,Wi)+h/2*F(T(i+1), Wi + hF(Ti,Wi))')
fprintf('\n');
d=input(' - Introduzca la ecuación diferencial encerrado con tildes y Dy sin xtiplicar a 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,'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 ',ya);
fprintf('\n');
for j=a:h:(b-h)
i=1+i;
w(i+1)=w(i)+((h/2)*(subs(f,{x,y},{t(i),w(i)})))+((h/2)*(subs(f,{x,y},{(t(i)+h),(w(i)+(h*(subs(f,{x,y},{t(i),w(i)}))))})));
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + h/2 f(t%1.0f,w%1.0f) + h/2 f(t%1.0f + h,w%1.0f + 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.9f,w%1.0f) + %1.5f f(%1.9f + %1.5f,w%1.0f + %1.5f f(%1.9f,w%1.0f))',i,i-1,h/2,t(i),i-1,h/2,t(i),h,i-1,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('\n');
end
fprintf('\n');
-------------------------------------------------------------------------------------------------------------
%MÉTODO DE EULER
fprintf('\n');
clear all
clc
fprintf(' ---------------\n')
fprintf(' MÉTODO DE EULER\n')
fprintf(' ---------------\n')
fprintf('\n');
syms x y
disp('Wi+1=Wi+hF(Ti,Wi)')
d=input(' - Introduzca la ecuación diferencial Dy no debe contener a y Dy+x=y y encerrarlo con : ');
n=input(' - Introduzca la condición y(a)=b encerrado con '' : ');
l=input(' - Introduzca la variable enerrado con '' : ');
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,l);
pretty(m);
fprintf('\n\n\n\n');
fprintf(' i ti wi y(t)');
%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
%Este for obtiene los valores aproximados de solución
for j=a:h:(b-h)
i=1+i;
w(i+1)=w(i)+(h*(subs(f,{x,y},{t(i),w(i)})));
end
%Presentación de los datos en forma de tabla utilizando matrices
[k' t' w' v']
-------------------------------------------------------------------------------------
%metodo de euler pero es archivo .m
function EULER(a,b,h,w0)
syms y t
ti=a;
w(1)=w0;
sol_exacta=dsolve('Dy=(t+y-1)^2','y(0)=2','t');
f=(t+y-1)^2;
n=(b-a)/h;
fprintf(' ti= %3.5f\n W=%5.10f \n\n',a,w0);
for i=0:n-1,
w(i+2)=w(i+1)+(h)*subs(f,{t,y},{ti,w(i+1)});
fprintf(' ti= %3.5f\n W= %5.10f \n\n',ti+h,w(i+2));
ti=ti+h;
end
syms ti wi
for i=0:0,
wi(i+2)=wi(i+1)+(h)*subs(f,{t,y},{ti,wi(i+1)});
pretty(wi(i+2));
end
pretty(sol_exacta);
ti=a;
for i=0:n,
z=subs(sol_exacta,t,ti);
fprintf('t=%3.5f Y= %5.10f \n',ti,z);
ti=ti+h;
end
------------------------------------------------------------------
metodo RK 4 para sistemas
disp('Método de solucion de sistemas RK 4 orden para 3 ecuaciones diferenciales');
m=input('Ingrese el # de ecuaciones: ');
%agregar mas variables simbolicas si hay mas ecuaciones
syms x y z t
fi=cell(1,m);
for i=1:m
fprintf('Ingrese f%d´: ',i); % 1+1/2*exp(t)-2*x
fi{i}=input('');
end
%fprintf('\n');
wi0=zeros(1,m);
for i=1:m
fprintf('Ingrese w%d0: ',i); % 2
wi0(i)=input('');
end
%fprintf('\n');
interv=input('Ingrese el intervalo: '); %[0 0.3]
h=input('Ingrese el tamaño de paso: '); % 0.1
n=(interv(2)-interv(1))/h;
n=round(n);
ti=zeros(1,n+1);
for i=2:n
ti(i)=interv(1)+(i-1)*h;
end
ti(n+1)=interv(2);
k1i=zeros(1,m);
k2i=zeros(1,m);
k3i=zeros(1,m);
k4i=zeros(1,m);
wij=zeros(n);
fprintf('\n');
for j=1:n
fprintf('\nj = %d\nt%d = %d\n\n',j-1,j-1, ti(j+1));
%fprintf('\nj=%d\n\n',j-1);
for i=1:m
if j==1
k1i(i)=h*subs(fi{i},{t,x,y,z},{ti(j),wi0(1),wi0(2),wi0(3)});
else
k1i(i)=h*subs(fi{i},{t,x,y,z},{ti(j),wij(1,j-1),wij(2,j-1),wij(3,j-1)});
end
%fprintf('K1%d = h*f(t%d,w%d)',i,j,j); %rk4 normalito
fprintf('k1%d = h*f%d(t%d',i,i,j-1);
for cont=1:m
fprintf(',w%d%d',cont,j-1);
end
fprintf(') = %.9f\n',k1i(i));
%fprintf('K1%d=%.9f\n',i,k1i(i)); %como lo tenia el gordo
end
for i=1:m
if j==1
k2i(i)=h*subs(fi{i},{t,x,y,z},{ti(j)+h/2,wi0(1)+1/2*k1i(1),wi0(2)+1/2*k1i(2),wi0(3)+1/2*k1i(3)});
else
k2i(i)=h*subs(fi{i},{t,x,y,z},{ti(j)+h/2,wij(1,j-1)+1/2*k1i(1),wij(2,j-1)+1/2*k1i(2),wij(3,j-1)+1/2*k1i(3)});
end
%fprintf('K2%d = h*f(t%d+h/2,w%d+1/2*k1)',i,j,j); %rk4 normalito
fprintf('k2%d = h*f%d(t%d+h/2',i,i,j-1);
for cont=1:m
fprintf(',w%d%d+(1/2)*k1%d',cont,j-1,cont);
end
fprintf(') = %.9f\n',k2i(i));
%fprintf('K2%d=%.9f\n',i,k2i(i));
end
for i=1:m
if j==1
k3i(i)=h*subs(fi{i},{t,x,y,z},{ti(j)+h/2,wi0(1)+1/2*k2i(1),wi0(2)+1/2*k2i(2),wi0(3)+1/2*k2i(3)});
else
k3i(i)=h*subs(fi{i},{t,x,y,z},{ti(j)+h/2,wij(1,j-1)+1/2*k2i(1),wij(2,j-1)+1/2*k2i(2),wij(3,j-1)+1/2*k2i(3)});
end
%fprintf('K3%d = h*f(t%d+h/2,w%d+1/2*k2)',i,j,j); %rk4 normalito
fprintf('k3%d = h*f%d(t%d+h/2',i,i,j-1);
for cont=1:m
fprintf(',w%d%d+(1/2)*k2%d',cont,j-1,cont);
end
fprintf(') = %.9f\n',k3i(i));
%fprintf('K3%d=%.9f\n',i,k3i(i));
end
for i=1:m
if j==1
k4i(i)=h*subs(fi{i},{t,x,y,z},{ti(j)+h,wi0(1)+k3i(1),wi0(2)+k3i(2),wi0(3)+k3i(3)});
else
k4i(i)=h*subs(fi{i},{t,x,y,z},{ti(j)+h,wij(1,j-1)+k3i(1),wij(2,j-1)+k3i(2),wij(3,j-1)+k3i(3)});
end
%fprintf('K4%d = h*f(t%d,w%d+k3)',i,j+1,j); %rk4 normalito
fprintf('k4%d = h*f%d(t%d+h',i,i,j-1);
for cont=1:m
fprintf(',w%d%d+k3%d',cont,j-1,cont);
end
fprintf(') = %.9f\n',k4i(i));
%fprintf('K4%d=%.9f\n',i,k4i(i));
end
for i=1:m
if j==1
wij(i,j)=wi0(i)+1/6*(k1i(i)+2*k2i(i)+2*k3i(i)+k4i(i));
else
wij(i,j)=wij(i,j-1)+1/6*(k1i(i)+2*k2i(i)+2*k3i(i)+k4i(i));
end
fprintf('w%d%d = w%d%d+(k1%d+2*k2%d+2*k3%d+k4%d)/6 = %.9f\n',i,j,i,j-1, i,i,i,i, wij(i,j));
%fprintf('w%d%d=%.9f\n',i,j,wij(i,j));
end
fprintf('\n');
-----------------------------------------------------------------------------------------------
%MÉTODO EXPLÍCITO DE ADAMS-BASHFORTH DE CUATRO PASOS
% - Introduzca la ecuación diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b : 'y(0)=0.5'
% - Introduzca la función de trabajo : y-(x^2)+1
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tamaño de paso h : 0.1
% - Introduzca la condición inicial : 0.5
fprintf('\n');
clear all
clc
fprintf(' ------------------------------------------------------\n')
fprintf(' MÉTODO DE EXPLÍCITO DE ADAMS BASHFORTH DE CUATRO PASOS\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 : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tamaño de paso h : ');
ya=input(' - Introduzca la condición inicial : ');
tx(1)=a;
vx(1)=a;
ax=a;
bx=b;
hx=h;
dx=0;
mx = dsolve(d,n,'x');
for px=ax:hx:bx
dx=1+dx;
tx(dx)=px;
vx(dx)=subs(mx,px);
end
yb=subs(mx,tx(2));
yc=subs(mx,tx(3));
yd=subs(mx,tx(4));
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;
w(2)=yb;
w(3)=yc;
w(4)=yd;
i=3;
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
a=a+(4*h);
k3=k(end);
fprintf('------------------------------------------------------------------------------------------------------------------');
fprintf('\n');
fprintf(' FÓRMULAS DE CADA ITERACIÓN');
fprintf('\n');
fprintf('------------------------------------------------------------------------------------------------------------------');
fprintf('\n\n');
fprintf('- w0 = %1.9f ',ya);
fprintf('\n\n');
fprintf('- w1 = %1.9f ',yb);
fprintf('\n\n');
fprintf('- w2 = %1.9f ',yc);
fprintf('\n\n');
fprintf('- w3 = %1.9f ',yd);
fprintf('\n');
%Este for obtiene los valores aproximados de solución
for j=a:h:b
i=1+i;
w(i+1)=w(i)+((h/24)*((55*(subs(f,{x,y},{t(i),w(i)})))-(59*(subs(f,{x,y},{t(i-1),w(i-1)})))+(37*(subs(f,{x,y},{t(i-2),w(i-2)})))-(9*(subs(f,{x,y},{t(i-3),w(i-3)})))));
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + h/24 [55*f(t%1.0f,w%1.0f) - 59*f(t%1.0f,w%1.0f) + 37*f(t%1.0f,w%1.0f) - 9*f(t%1.0f,w%1.0f)]',i,i-1,i-1,i-1,i-2,i-2,i-3,i-3,i-4,i-4);
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + %1.9f [55*f(%1.9f,w%1.0f) - 59*f(%1.9f,w%1.0f) + 37*f(%1.9f,w%1.0f) - 9*f(%1.9f,w%1.0f)]',i,i-1,h/12,t(i),i-1,t(i-1),i-2,t(i-2),i-3,t(i-3),i-4);
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.9f %10.9f %10.9f',k(k2),t(k2),w(k2),v(k2));
fprintf('\n');
end
fprintf('\n');
-----------------------------------------------------------------------------------------------
%MÉTODO EXPLÍCITO DE ADAMS-BASHFORTH DE CINCO PASOS
% - Introduzca la ecuación diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condición y(a)=b : 'y(0)=0.5'
% - Introduzca la función de trabajo : y-(x^2)+1 despeje Dy
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tamaño de paso h : 0.1
% - Introduzca la condición inicial : 0.5
% x=w1j y=w2j
fprintf('\n');
clc
fprintf(' -----------------------------------------------------\n')
fprintf(' MÉTODO DE EXPLÍCITO DE ADAMS BASHFORTH DE CINCO PASOS\n')
fprintf(' -----------------------------------------------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuación diferencial (x,y) : ');
n=input(' - Introduzca la condición y(a)=b : ');
f1=input(' - Introduzca la función de trabajo (x,y) : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tamaño de paso h : ');
ya=input(' - Introduzca la condición inicial : ');
tx(1)=a;
vx(1)=a;
ax=a;
bx=b;
hx=h;
dx=0;
mx = dsolve(d,n,'x');
for px=ax:hx:bx
dx=1+dx;
tx(dx)=px;
vx(dx)=subs(mx,px);
end
yb=subs(mx,tx(2));
yc=subs(mx,tx(3));
yd=subs(mx,tx(4));
ye=subs(mx,tx(5));
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;
w(2)=yb;
w(3)=yc;
w(4)=yd;
w(5)=ye;
i=4;
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
a=a+(5*h);
k3=k(end);
fprintf('---------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n');
fprintf(' FÓRMULAS DE CADA ITERACIÓN');
fprintf('\n');
fprintf('---------------------------------------------------------------------------------------------------------------------------------------------------');
fprintf('\n\n');
fprintf('- w0 = %1.9f ',ya);
fprintf('\n\n');
fprintf('- w1 = y(t1)= %1.15f ',yb);
fprintf('\n\n');
fprintf('- w2 = y(t2)= %1.15f ',yc);
fprintf('\n\n');
fprintf('- w3 = y(t3)=%1.15f ',yd);
fprintf('\n\n');
fprintf('- w4 = y(t4)=%1.15f ',ye);
fprintf('\n\n');
%Este for obtiene los valores aproximados de solución
for j=a:h:b
i=1+i;
w(i+1)=w(i)+((h/720)*((1901*(subs(f,{x,y},{t(i),w(i)})))-(2774*(subs(f,{x,y},{t(i-1),w(i-1)})))+(2616*(subs(f,{x,y},{t(i-2),w(i-2)})))-(1274*(subs(f,{x,y},{t(i-3),w(i-3)})))+(251*(subs(f,{x,y},{t(i-4),w(i-4)})))));
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + h/720 [1901*f(t%1.0f,w%1.0f) - 2774*f(t%1.0f,w%1.0f) + 2616*f(t%1.0f,w%1.0f) - 1274*f(t%1.0f,w%1.0f) + 251*f(t%1.0f,w%1.0f)]',i,i-1,i-1,i-1,i-2,i-2,i-3,i-3,i-4,i-4,i-5,i-5);
fprintf('\n');
fprintf('- w%1.0f = w%1.0f + %1.9f [1901*f(%1.9f,w%1.0f) - 2774*f(%1.9f,w%1.0f) + 2616*f(%1.9f,w%1.0f) - 1274*f(%1.9f,w%1.0f) + 251*f(%1.9f,w%1.0f)]',i,i-1,h/12,t(i),i-1,t(i-1),i-2,t(i-2),i-3,t(i-3),i-4,t(i-4),i-5);
fprintf('\n');
fprintf(' w%1.0f = %4.15f ',i,w(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('\n');
end
fprintf('\n');
Suscribirse a:
Entradas (Atom)