Trying to get my runge-kutta to plot two graphs one of the original differential and one exact solution
8 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Abraham Boayue
2022-2-13
编辑:Abraham Boayue
2022-2-13
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 个评论
Abraham Boayue
2022-2-13
The the integration was done correctly. Your asnswer for the exact solution is right.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!