Trying to get my runge-kutta to plot two graphs one of the original differential and one exact solution
16 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Zachary David
il 13 Feb 2022
Commentato: Abraham Boayue
il 13 Feb 2022
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
0 Commenti
Risposta accettata
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
Abraham Boayue
il 13 Feb 2022
The the integration was done correctly. Your asnswer for the exact solution is right.
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!