%_________________________________________________________________________Ode para aproximada (Cualquier ode)
%[t,y]=ode(function,[a:h:b],cond_inicial)
%fdy=archivo m, y(pi/4)=0; pi/4<=t<=5pi*12; h=pi/50
format long
[t,y]=ode23(@fdy,[pi/4:pi/50:5*pi/12],0)%Valor aproximado
ED=\'Dy=2*sin(2*t)-1-2*y\';
cond_inicial=\'y(pi/4)=0\';
exacta=dsolve(ED,cond_inicial) %Valor exacto
Y=double(subs(exacta,t))
error=norm(Y-y,inf)
El archivo M es:
function dy=fdy(t,y)
dy=2*sin(2*t)-1-2*y;
end
%____________________________________________________________________________ Sistema primer orden
%Archivo m
function dy=fdy(t,xy)
%escribir xy en vez de solo x o y, pero poner (1) para X and (2) para Y
dy=[(1/8)*(exp(t)+1+t-17*xy(1));(1/8)*(8*exp(t)-3-6*t-17*xy(2))];
end
%[t xy]=ode45(@fdy,[a:h:b],[cond_inicial X;cond_inicial Y]);
format long
[t xy]=ode45(@fdy,[0:0.1:1],[1;0])
%Hallando valor exacto
EDX=\'Dx+17*x/8=(1/8)*(exp(t)+1+t)\';%Cambiar acá
EDY=\'Dy+17*y/8=(1/8)*(8*exp(t)-3-6*t)\';%Cambiar acá
xi=\'x(0)=1\';%Cambiar acá
yi=\'y(0)=0\';%cambiar acá
[x y]=dsolve(EDX,EDY,xi,yi)
exacta_x=double(subs(x,t));
exacta_y=double(subs(y,t));
exacta=[exacta_x exacta_y]
errorx=norm(exacta_x-xy(:,1));
errory=norm(exacta_y-xy(:,2));
error=[errorx errory]
%______________________________________________________________________________Orden Superior
%Archivo m
function du=fdy(t,u)
du=zeros(3,1);%Cambiar el 3 por 4 si es de orden 4
du(1)=u(2);
du(2)=u(3);
du(3)=exp(t)-2*u(3)+u(2)+2*u(1); (La ecuación de acá se obtiene al despejar y^n)
end
%Si fuera de orden 4, solo agregar du(4)=D4y
%Para ese caso sería du(3)=u(4)
format long
%ode(@f,[a:h:b],[ciY;ciDy;ciD2y])
[t,XD]=ode45(@fdy,[0:0.05:0.5],[1;2;0])%Valor aproximado
%valor exacto
exacta=dsolve(\'D3y+2*D2y-Dy-2*y=exp(t)\',\'y(0)=1\',\'Dy(0)=2\',\'D2y(0)=0\')
Y=double(subs(exacta,t))
%error
error=norm(Y-XD(:,1),inf)
%______________________________________________________________________________ Euller modificado
function M=EulerMod(a,b,h,ci,f,fd,cied)
syms t y
format long
n=int32((b-a)/h);
fexacta = dsolve(fd,cied);
ti=a:h:b;
wi(1)=ci;
exacta=double(subs(fexacta,t,ti));
for i=1:n
wi(i+1)=wi(i)+h/2*(subs(f,[t,y],[ti(i),wi(i)]))+h/2*(subs(f,[t,y],[ti(i+1),wi(i)+h*(subs(f,[t,y],[ti(i),wi(i)]))]));
end
M=[ti\' wi\' exacta\'];
error=norm(M(:,3)-M(:,2),inf)
fprintf(\'\\tT \\t\\tWi\\t\\tY(t)\')
%f=función de trabajo (\'t^2+5y\')
%cied=y(k)=alpha; \'y(0)=1\';
%fd=función pero con el diferencial (\'Dy=t^2+5y\')
%ci=cond_inicial (alpha=1)
%___________________________________________________________________________________Heun
function M=Heun(a,b,h,ci,f,ed,ya)
% Heun(0,1/50,1/200,1/3,\'5*t^2+2*t-5*y\',\'Dy=5*t^2+2*t-5*y\',\'y(0)=1/3\')
syms t y
n=int32((b-a)/(h));
ti=a:h:b;
fexacta = dsolve(ed,ya);
exacta=double(subs(fexacta,t,ti));
wi(1)=ci;
for i=1:n
wi(i+1)=wi(i)+h/4*(subs(f,[t,y],[ti(i),wi(i)]))+(3*h)/4*(subs(f,[t,y],[ti(i)+(2*h)/3,wi(i)+(2*h)/3*(subs(f,[t,y],[ti(i),wi(i)]))]));
end
fprintf(\' ti wi exacta \');
fprintf(\'\\n\');
M=[ti\' wi\' exacta\'];
%_________________________________________________________________________________ Euller
function M=Euler(a,b,h,ci,f,fd,cied)
%Euler(1,61/60,1/300,10,\'(1/(t+1))*(log(t)-y)\',\'Dy=(1/(t+1))*(log(t)-y)\',\'y(1)=10\')
syms t y
n = int32(((b-a)/h));
fexacta = dsolve(fd,cied);
ti=a:h:b;
wi(1)=ci;
exacta = double(subs(fexacta,t,ti));
for i=1:n
wi(i+1)=wi(i)+h*(subs(f,[t,y],[ti(i),wi(i)]));
end
M=[ti\' wi\' exacta\'];
%_____________________________________________________________________________________Punto medio
function M=PuntoMedio(a,b,h,ci,f,ed,ya)
%PuntoMedio(1,9/8,1/80,2,\'(y^3-t^3)/(t*y^2)\',\'Dy=(y^3-t^3)/(t*y^2)\',\'y(1)=2\')
syms t y
n=int32((b-a)/h);
ti=a:h:b;
fexacta = dsolve(ed,ya);
exacta=double(subs(fexacta,t,ti));
wi(1)=ci;
for i=1:n
wi(i+1)=wi(i)+h*(subs(f,[t,y],[ti(i)+h/2,wi(i)+h/2*(subs(f,[t,y],[ti(i),wi(i)]))]));
end
fprintf(\' ti wi exacta \');
fprintf(\'\\n\');
M=[ti\' wi\' exacta\'];
%______________________________________________________________________________________Runge-Kutta-Orden4
function M=RKOrden4(a,b,h,ci,f,fd,cied)%RUNGE-KUTTA ORDEN 4
syms t y
format long
n=int32((b-a)/h);
fexacta = dsolve(fd,cied);
ti=a:h:b;
wi(1)=ci;
exacta=double(subs(fexacta,t,ti));
for i=1:n
k1=h*(subs(f,[t,y],[ti(i),wi(i)]));
k2=h*(subs(f,[t,y],[ti(i)+h/2,wi(i)+1/2*k1]));
k3=h*(subs(f,[t,y],[ti(i)+h/2,wi(i)+1/2*k2]));
k4=h*(subs(f,[t,y],[ti(i+1),wi(i)+k3]));
wi(i+1)=wi(i)+1/6*(k1+2*k2+2*k3+k4);
end
M=[ti\' wi\' exacta\'];
%Caso de aplicación del circuito RLC
%Ly\'\' + Ry\' +y/c = E ___________ q(a)=k; i(a)=k