Amor mío ED

%_________________________________________________________________________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 ttWittY(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