Trying to get my runge-kutta to plot two graphs one of the original differential and one exact solution

16 visualizzazioni (ultimi 30 giorni)
This is what I have so far but it just sits there and doesnt graph anything. Is there something else I need to do or is there something I am missing?
clear,clc
function [x,y] = rk4(ti,tf,npts,yi,xi,F_tx,F_ty)
h=0.1;
ti = 0;
tf = 24;
t=ti:h:tf;
xi= 50;
yi= 50;
F_tx=@(x,y)(75-75*exp^(-.25*t));
F_ty=@(x,y)(0.25*(75-yi));
x=zeros(1,length(t));
x(1)=xi;
y=zeros(1,length(t));
y(1)=yi;
%t= linspace(ti,tf,length(t));
for i=1:(length(t)-1)
kx1 = F_tx(x(i),y(i));
ky1 = F_ty(x(i),y(i));
kx2 = F_tx(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
ky2 = F_ty(x(i)+0.5*h*kx1,y(i)+0.5*h*ky1);
kx3 = F_tx((x(i)+0.5*h*kx2),y(i)+0.5*h*ky2);
ky3 = F_ty((x(i)+0.5*h*kx2),(y(i)+0.5*h*ky2));
kx4 = F_tx((x(i)+kx3*h),y(i)+ky3*h);
ky4 = F_ty((x(i)+kx3*h),(y(i)+ky3*h));
x(i+1) = x(i) + h*(kx1+2*kx2+2*kx3+kx4)/6;
y(i+1) = y(i) + h*(ky1+2*ky2+2*ky3+ky4)/6;
end
figure(1)
plot(t,y)
hold on
plot(t,x)
end

Risposta accettata

Abraham Boayue
Abraham Boayue il 13 Feb 2022
Modificato: Abraham Boayue il 13 Feb 2022
So you meant to solve dT/dt = 0.25(75-T) with initial condition of T(0) = 0? If so, then the exaction solution is T(t) = 75-75e^(-0.25t). Let me know if this is what you wanted.
function [y, t,N] = rk4(f,to,tf,ya,h)
N= ceil((tf-to)/h);
halfh = h / 2;
y(1,:) = ya;
t(1) = to;
h6 = h/6;
for i = 1 : N
t(i+1) = t(i) + h;
th2 = t(i) + halfh;
s1 = f(t(i), y(i,:));
s2 = f(th2, y(i,:) + halfh * s1);
s3 = f(th2, y(i,:) + halfh * s2);
s4 = f(t(i+1), y(i,:) + h * s3);
y(i+1,:) = y(i,:) + (s1 + s2+s2 + s3+s3 + s4) * h6;
end
%% This is your mfile. Run this alone after you have saved the above function in the same
% folder as the mfile.
yo = 0;
h = 0.1;
to = 0;
tf = 24;
Exact =@(t)(75-75*exp(-.25*t));
f = @(t,y)(0.25*(75-y));
[y, t,N] = rk4(f,to,tf,yo,h);
plot(t,Exact(t),'linewidth',2.5)
hold on
plot(t,y,'--','linewidth',2.5)
b = xlabel('t');
set(b,'Fontsize',15);
b = ylabel('y');
set(b,'Fontsize',15);
a= title('Exact Solution and Numerical solution');
set(b,'Fontsize',15);
legend('Exact Soln.','Numerical Soln.','location','best')
grid
  3 Commenti

Accedi per commentare.

Più risposte (0)

Prodotti


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by