How do I fix this runge-kutta algorithm implimentation?
显示 更早的评论
I have a matrix c with 3 columns containing 3 sets of data I want to plot, mimicing 4 bolus injections total, changing the "q" value to modify blood concentration at each injection to maintain a certain concentration level. I used the following documentation as a guide to establish the algorithm itself but I have reached a standstill: https://www.mathworks.com/videos/solving-odes-in-matlab-3-classical-runge-kutta-ode4-117528.html
%4 bolues total - at time 0, 24, 48, 72
%1st function - dcdt=f-mc
%this is dydt in the documentation
dcdt=@(t,c,f)(f-Mc);
%impliment algorithm
%pick parameters of start and stop time and step size
%solve for what f should be in each step of time
%use create_f, give it t, h, q, feed it into dcdt to calculate k values
%update c and t
t0=0; %start time
tfinal=94; %stopped time
h=0.1; %step size
%initial conditions
c=1;
c0=0;
t=(0:h:tfinal)';
%dcdt=create_f(t,h,5)
something=ode4(dcdt,t0,h,tfinal,c0);
plot(t,something)
%runge-kutta algorithm implimentation
function cout=ode4(dcdt,t0,h,tfinal,c0)
c=c0;
cout=c;
for t=t0:h:tfinal-h
k1=dcdt(t, c);
k2=dcdt(t+(h/2), c+h*(k1/2));
k3=dcdt(t+(h/2), c+h*(k2/2));
k4=dcdt(t+h, c+h*k3);
c=c+h*(k1+2*k2+2*k3+k4)/6;
cout=[cout; c];
end
end
%2nd function - define function to create f based on time
%f = [f1; f2; f3], f2=f3=0 for all time
%f1=sum of boluses
%must take time, step size, q which is delta(t-tb), q is plasma []
%immediately after bolus injuction. play with q until i get the target concentrations within the bounds
function out=create_f(t,h,q)
f1=(q/h)*((u(t)-u(t-h))+(u(t-24)-u(t-24-h))+(u(t-48)-u(t-48-h))+(u(t-72)-u(t-72-h)));
f2=0;
f3=0;
out=[f1; f2; f3];
end
%3rd function - unit step, takes time, returns output of unit step function out=u(t)
if t>=0
out=1;
else
out=0;
end
end
2 个评论
Jan
2023-2-6
Please explain, what "I have reached a standstill" means. You do know, what the problem is, but the readers don't. Instead of spending time to help you, they have to guess at first, what you consider as problem.
Torsten
2023-2-6
dcdt=@(t,c,f)(f-Mc);
M is not defined.
And most probably you mean M*c instead of Mc.
k1=dcdt(t, c);
k2=dcdt(t+(h/2), c+h*(k1/2));
k3=dcdt(t+(h/2), c+h*(k2/2));
k4=dcdt(t+h, c+h*k3);
It will work, but you should notice that dcdt is defined with three arguments (t,c,f).
回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Biological and Health Sciences 的更多信息
2023-2-6
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
