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
Best Answer