algunos metodos

Esta es una vista previa del archivo. Inicie sesión para ver el archivo original

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蒚ODO 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
Publicado por M閠odos en matlab en 15:13 18 comentarios: 
Enviar por correo electr髇ico
Escribe un blog
Compartir con Twitter
Compartir con Facebook
Compartir en Pinterest
-----------------------------------------------------------------------------------------
%M蒚ODO DEL PUNTO MEDIO
fprintf('\n');
clear all
clc
fprintf(' ----------------------\n')
fprintf(' M蒚ODO 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髇 diferencial encerradas en tildes y q Dy no este x y : ');
n=input(' - Introduzca la condici髇 y(a)=b encerradas en tildes : ');
l=input(' - Introduzca la variable encerradas en tildes : ');
f1=input(' - Introduzca la funci髇 de trabajo : ');
ya=input(' - Introduzca la condici髇 inicial : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tama駉 de paso h : ');
fprintf('\n\n');
fprintf(' - La soluci髇 de la ecuaci髇 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閚 se utiliza para evaluar la ecuaci髇 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髇
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髇 de los datos en forma de tabla utilizando matrices
[k' t' w' v'] 
----------------------------------------------------------------------------------------------------------------
%M蒚ODO DE HEUN
fprintf('\n');
clear all
clc
fprintf(' --------------\n')
fprintf(' M蒚ODO 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髇 diferencial encerrados en tildes y q DY no est multiplpcd * y : ');
n=input(' - Introduzca la condici髇 y(a)=b : ');
l=input(' - Introduzca la variable : ');
f1=input(' - Introduzca la funci髇 de trabajo : ');
ya=input(' - Introduzca la condici髇 inicial : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tama駉 de paso h : ');
fprintf('\n\n');
fprintf(' - La soluci髇 de la ecuaci髇 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閚 se utiliza para evaluar la ecuaci髇 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髇
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髇 de los datos en forma de tabla utilizando matrices
[k' t' w' v'] 
_-----------------------------------------------------------------------------------------------------------
%M蒚ODO DE RUNGE-KUTTA-FEHLBERG ORDEN 5
% - Introduzca la ecuaci髇 diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condici髇 y(a)=b : 'y(0)=0.5'
% - Introduzca la funci髇 de trabajo : y-(x^2)+1
% - Introduzca la condici髇 inicial : 0.5
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tama駉 de paso h : 0.1
clear all
disp('---------------------------------------------------------------------------');
disp('M閠odo de RUNGE-KUTTA-FEHLBERG de orden 5');
disp('---------------------------------------------------------------------------');
fprintf('\n')
syms x y
d=input(' - Ingrese la ecuaci髇 diferencial en (x,y): ');
n=input(' - Ingrese la condici髇 y(a)=b: ');
f1=input(' - Ingrese la funci髇 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髇 de la ecuaci髇 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髇: %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閚 se utiliza para evaluar la ecuaci髇 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髇 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駉 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蒚ODO 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髇 diferencial encerrado con tildes y Dy sin xtiplicar a y : ');
n=input(' - Introduzca la condici髇 y(a)=b : ');
f1=input(' - Introduzca la funci髇 de trabajo : ');
ya=input(' - Introduzca la condici髇 inicial : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tama駉 de paso h : ');
fprintf('\n\n');
fprintf(' - La soluci髇 de la ecuaci髇 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閚 se utiliza para evaluar la ecuaci髇 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髇
fprintf('-------------------------------------------------------------------------------------------------------');
fprintf('\n');
fprintf(' F覴MULAS DE CADA ITERACI覰');
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髇 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蒚ODO DE EULER
fprintf('\n');
clear all
clc
fprintf(' ---------------\n')
fprintf(' M蒚ODO DE EULER\n')
fprintf(' ---------------\n')
fprintf('\n');
syms x y
disp('Wi+1=Wi+hF(Ti,Wi)')
d=input(' - Introduzca la ecuaci髇 diferencial Dy no debe contener a y Dy+x=y y encerrarlo con : ');
n=input(' - Introduzca la condici髇 y(a)=b encerrado con '' : ');
l=input(' - Introduzca la variable enerrado con '' : ');
f1=input(' - Introduzca la funci髇 de trabajo : ');
ya=input(' - Introduzca la condici髇 inicial : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tama駉 de paso h : ');
fprintf('\n\n');
fprintf(' - La soluci髇 de la ecuaci髇 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閚 se utiliza para evaluar la ecuaci髇 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髇
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髇 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閠odo 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駉 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蒚ODO EXPL虲ITO DE ADAMS-BASHFORTH DE CUATRO PASOS
% - Introduzca la ecuaci髇 diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condici髇 y(a)=b : 'y(0)=0.5'
% - Introduzca la funci髇 de trabajo : y-(x^2)+1
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tama駉 de paso h : 0.1
% - Introduzca la condici髇 inicial : 0.5
fprintf('\n');
clear all
clc
fprintf(' ------------------------------------------------------\n')
fprintf(' M蒚ODO DE EXPL虲ITO DE ADAMS BASHFORTH DE CUATRO PASOS\n')
fprintf(' ------------------------------------------------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuaci髇 diferencial : ');
n=input(' - Introduzca la condici髇 y(a)=b : ');
f1=input(' - Introduzca la funci髇 de trabajo : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tama駉 de paso h : ');
ya=input(' - Introduzca la condici髇 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髇 de la ecuaci髇 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閚 se utiliza para evaluar la ecuaci髇 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覴MULAS DE CADA ITERACI覰');
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髇
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髇 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蒚ODO EXPL虲ITO DE ADAMS-BASHFORTH DE CINCO PASOS
% - Introduzca la ecuaci髇 diferencial : 'Dy=y-(x^2)+1'
% - Introduzca la condici髇 y(a)=b : 'y(0)=0.5'
% - Introduzca la funci髇 de trabajo : y-(x^2)+1 despeje Dy
% - Introduzca el valor de a : 0
% - Introduzca el valor de b : 1
% - Introduzca el tama駉 de paso h : 0.1
% - Introduzca la condici髇 inicial : 0.5
% x=w1j y=w2j
fprintf('\n');
clc
fprintf(' -----------------------------------------------------\n')
fprintf(' M蒚ODO DE EXPL虲ITO DE ADAMS BASHFORTH DE CINCO PASOS\n')
fprintf(' -----------------------------------------------------\n')
fprintf('\n');
syms x y
d=input(' - Introduzca la ecuaci髇 diferencial (x,y) : ');
n=input(' - Introduzca la condici髇 y(a)=b : ');
f1=input(' - Introduzca la funci髇 de trabajo (x,y) : ');
a=input(' - Introduzca el valor de a : ');
b=input(' - Introduzca el valor de b : ');
h=input(' - Introduzca el tama駉 de paso h : ');
ya=input(' - Introduzca la condici髇 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髇 de la ecuaci髇 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閚 se utiliza para evaluar la ecuaci髇 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覴MULAS DE CADA ITERACI覰');
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髇
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髇 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');
 
Publicado por M閠odos en matlab en 14:42 No hay comentarios: 
Enviar por correo electr髇ico
Escribe un blog
Compartir con Twitter
Compartir con Facebook
Compartir en Pinterest
jueves, 3 de abril de 2014
% integraci髇 de romber
clear all;
syms x y z q valorsitos;
disp('Metodo de Romberg')
f=input('Ingrese la funcion: ');
a=input('Ingrese el valor de a: ');
b=input('Ingrese el valor de b: ');
n=input('Ingrese el valor de N: ');
cont1=1;
%encontrando los valores de h
while(cont1<=n)
 h(cont1)=(b-a)/(2^(cont1-1));
 fprintf('h%d=%.9f\n\n',cont1,h(cont1));
 cont1=cont1+1;
end
%generando la primera columna
t(1,1)=(h(1)/2)*(subs(f,a)+subs(f,b));
fprintf('R(1,1)=h1/2*(f(a)+f(b))==>');
disp(t(1,1));
cont1=2;
sumatoria=0;
while(cont1<=n)
 cont2=1;
 fprintf('R(%1.0f,%1.0f)=',cont1,1);
 valorsitos(cont1,cont2)=q+h(cont1-1);
 while(cont2<=(2^(cont1-2)))
 
 valorsitos(cont1,cont2+1)=(q+a+(2*cont2-1)*h(cont1));
 sumatoria=sumatoria+subs(f,(a+(2*cont2-1)*h(cont1)));
 cont2=cont2+1;
 end
 t(cont1,1)=(1/2)*(t(cont1-1,1)+h(cont1-1)*sumatoria);
 fprintf(')==>');
 disp(t(cont1,1));
 sumatoria=0;
 cont1=cont1+1;
end
%generando el resto de la tabla
cont1=2;
cont2=2;
cont3=2;
while(cont2<=n)
 while(cont1<=n)
 t(cont1,cont2)=t(cont1,cont2-1)+(t(cont1,cont2-1)-t(cont1-1,cont2-1))/(4^(cont2-1)-1);
 fprintf('R(%1.0f,%1.0f)=R(%1.0f,%1.0f)+(R(%1.0f,%1.0f)-R(%1.0f,%1.0f)/%d)==>',cont1,cont2,cont1,cont2-1,cont1,cont2-1,cont1-1,cont2-1,(4^(cont2-1)-1));
 disp(t(cont1,cont2));
 cont1=cont1+1;
 end
 cont3=cont3+1;
 cont1=cont3;
 cont2=cont2+1;
end
%clear all integracion de simpsom
syms x y;
funcion=input('funcion=');
a=input('a=');
b=input('b=');
n=input('n=');
h=(b-a)/n;
for i=1:n-1
 x(i)=a+i*h;
 fprintf('x(%d)= ',i)
 pretty(x(i));
 %disp(x(i));
end
part2=0;
for i=1:(n/2-1)
 part2=part2+subs(funcion,x(i*2));
end
part3=0;
for i=1:(n/2)
 part3=part3+subs(funcion,x(i*2-1));
end
resultado=(h/3)*(subs(funcion,a)+2*part2+4*part3+subs(funcion,b));
double(resultado)
 % extrapolaci髇 de richarson programa 1
function N = richardson2B(f, x0, n, h, pre)
%meter lo de arriba asi como va copy paste! y la funcion f es primero
%hayq declararla!
% f es la funcion, x0 el de aproximar, n el subindice de N, h es h y pre es
% precisi锟絥 (opcional)
if nargin < 5
 pre = 8;
end
pre = num2str(pre);
N = zeros(n,n);
% generamos la primera columna
for a = 1:n
 h_ = h/(2^(a-1)); 
 expr = sprintf('[1/(2*(%g))]*[f(%g) - f(%g)]', h_, x0 + h_, x0 - h_); 
 N(a, 1) = centrada3_fun(f, x0, h_, str2num(pre)); 
 
end
% las demas columnas
for a = 2:n
 for b = a:n
 N(b, a) = N(b, a-1) + (N(b, a-1) - N(b-1, a-1))/(4^(a-1) - 1);
 end
end
 % extrpolacion de Richardson programa 2
function RICHARDSON(h,N,Xo,g)
%h es el espaciado, N es hasta q numero se va a llegar N3,N4, etc. Xo es el
%punto en donde se va a evaluar la derivada y g es la funcion.
R=zeros(N,N);
for j=1:N
 X1=Xo+h;
 X2=Xo-h;
 syms x
 g1=subs(g,X1);
 g2=subs(g,X2);
 %Primero usamos la formula de 3 puntos
 %f'(h)=(1/2h)(f(x0+h)-f(x0-h))
 %f'(h/2)=(1/2(h/2))(f(x0+(h/2))-f(x0-(h/2)))
 R(j,1)=(1/(2*h))*(g1-g2);
 fprintf('h=%1.4f\n',h)
 h=h/2;
end
for j=2:N
 for i=j:N
 R(i,j)=R(i,j-1) + (R(i,j-1)-R(i-1,j-1))/((4^(j-1))-1);
 end
end
for j=1:N
 for i=j:N
fprintf('N(%d,%d)= %3.10f\n',i,j,R(i,j))
 end
end
valor_aproximado=R(N,N)
gp=diff(g);
valor_exacto=subs(gp,Xo)
error=abs(valor_aproximado - valor_exacto)
%Formula de Richardson
% N1(h/2)-N1(h) 
%N2(h)=N1(h/2)+ -------------
% 4^(j-1))-1
%
% N1(h/4)-N1(h/2) 
%N2(h/2)=N1(h/4)+ -------------
% 4^(j-1))-1
%
%....