# How to make a function that uses Runge-Kutta Method

12 views (last 30 days)
H-M on 16 Jun 2020
Commented: darova on 27 Jun 2020
Here is my code to find the velocity and position of particles in 3 dimensions using Rung-Kutta 4th order method as a time marching:
I just mentined a part of code here,its a little long.
Uf=Poiseuille flow(fluid) in pipe and wp,up,vp are particles velocity in 3 dimensions and FW=dwp/dt,FU=dup/dt,FV=dvp/dt.
k1-k4 and m1-m4 and l1-l4 are slopes from Runge-Kutta 4th oreder method to obtain velocity and position of particles.
MY QUESTION is how can I read a Runge-Kutta 4th as an function to avoid repeating these k,m and l's and just call the that function.
while ( t(j)<Time )
FW=@(wp,Uf) (1/(rop+0.5*ro))*((rop-ro)*g+(18*miu/Dp^2)*(Uf-wp)); % wp is z's particle velocity
FU=@(up,Uf) (1/(rop+0.5*ro))*((rop-ro)*0+(18*miu/Dp^2)*(0 -up)); % up is x's particle velocity
FV=@(vp,Uf) (1/(rop+0.5*ro))*((rop-ro)*0+(18*miu/Dp^2)*(0 -vp)); % vp is y's particle velocity
k1 = FW(wp(:,:,:,j),Uf)*Dt; m1 = FU(up(:,:,:,j),Uf)*Dt; l1 = FV(vp(:,:,:,j),Uf)*Dt;
k2 = FW(wp(:,:,:,j)+0.5*k1,Uf)*Dt; m2 = FU(up(:,:,:,j)+0.5*m1,Uf)*Dt; l2 = FV(vp(:,:,:,j)+0.5*l1,Uf)*Dt;
k3 = FW(wp(:,:,:,j)+0.5*k2,Uf)*Dt; m3 = FU(up(:,:,:,j)+0.5*m2,Uf)*Dt; l3 = FV(vp(:,:,:,j)+0.5*l2,Uf)*Dt;
k4 = FW(wp(:,:,:,j)+k3,Uf)*Dt; m4 = FU(up(:,:,:,j)+m3,Uf)*Dt; l4 = FV(vp(:,:,:,j)+l3,Uf)*Dt;
wp(:,:,:,j+1) = wp(:,:,:,j) + (1/6)*(k1+2*k2+2*k3+k4); % particle velocity along the pipe.
up(:,:,:,j+1) = up(:,:,:,j) + (1/6)*(m1+2*m2+2*m3+m4);
vp(:,:,:,j+1) = vp(:,:,:,j) + (1/6)*(l1+2*l2+2*l3+l4);
Zp(:,:,:,j+1) = Zp(:,:,:,j)+Dt*wp(:,:,:,j); % particle displacement along the pipe.
Xp(:,:,:,j+1) = Xp(:,:,:,j)+Dt*up(:,:,:,j);
Yp(:,:,:,j+1) = Yp(:,:,:,j)+Dt*vp(:,:,:,j);
j=j+1;
t(j)=Dt*j;
end

darova on 17 Jun 2020
Edited: darova on 17 Jun 2020
Try something like that
wp(j+1) = rk4(FW,t,wp(j),h);
function u1 = rk4(f,t,u0,h)
k1 = f(t,u0);
k2 = f(t+h/2,u0+k1*h/2);
k3 = f(t+h/2,u0+k2*h/2);
k4 = f(t+h,u0+h*k3);
u1 = u0 + 1/6*h*(k1 + 2*k2 + 2*k3 + k4);
darova on 27 Jun 2020
i think something like that
for i = 1:length(t)
wp(j+1) = rk4(FW,t(i),wp(i),Dt);
end