Problem solving simultaneous ODEs using 4th order Runge-Kutta method (Too many input arguments.)
Show older comments
I am trying to solve 6 ODEs simulataneously (three 2nd order and three 1st order). I am solving it using 4th order Runge-Kutta method. But I am getting an error: Too many input arguments.
Please suggest. My code is following.
tic
clear;
close all;
clc;
ae=10 ; % aspect ratio
ds=1*10^-3 ; % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ; % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2; % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ; % time period
B=1.1;
theta_ini= 0 ; % azimuthal angle matrix
phi_ini= pi/4 ; % polar angle matrix
vp= 0.1*w1*A ; % Particle intrinsic velocity (m/s)
H=3700; % depth of ocean
Fr=w1*A./sqrt(g*H); % Froude number
z0= 3 ;
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*T; %seconds
dt=0.05 ;
t=t0:dt:tf ;
c1=18*nu/(B*ds^2) ;
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
% X=[x; vx; y; vy; z; vz; px; py; pz];
% RK4
fx = @(x,vx,y,vy,z,vz,px,py,pz) vx;
fvx = @(x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
(1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t))));
fy = @(x,vx,y,vy,z,vz,px,py,pz) vy;
fvy = @(x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*C22d)).*(-c1.*(K22.*vy));
fz = @(x,vx,y,vy,z,vz,px,py,pz) vz;
fvz = @(x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
(1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g);
fpx = @(x,vx,y,vy,z,vz,px,py,pz) ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpy = @(x,vx,y,vy,z,vz,px,py,pz) -ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpz = @(x,vx,y,vy,z,vz,px,py,pz) ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
h=dt;
x=zeros(length(t),1);
vx=zeros(length(t),1);
px=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
py=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
pz=zeros(length(t),1);
% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;
u0y=0;
u0z=0;
% Initial particle velocity
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;
% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);
x(1) = 0 ;
vx(1) = v_ini_1p ;
y(1) = 0 ;
vy(1) = v_ini_2p ;
z(1) = -z0 ;
vz(1) = v_ini_3p ;
px(1) = initial_p1 ;
py(1) = initial_p2 ;
pz(1) = initial_p3 ;
for i=1: length(t)-1
k_1x = fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1vx = fvx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1y = fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1vy = fvy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1z = fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1vz = fvz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
K1y = Fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
K1z = Fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1px = fpx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1py = fpy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1pz = fpz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_2x = fx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2y = fy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2z = fz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_3x = fx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3y = fy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3z = fz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_4x = fx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4y = fy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4z = fz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
x(i+1) = x(i) + (1/6)*(k_1x+2*k_2x+2*k_3x+k_4x)*h;
vx(i+1) = vx(i) + (1/6)*(k_1vx+2*k_2vx+2*k_3vx+k_4vx)*h;
px(i+1) = px(i) + (1/6)*(k_1px+2*k_2px+2*k_3px+k_4px)*h;
y(i+1) = y(i) + (1/6)*(k_1y+2*k_2y+2*k_3y+k_4y)*h;
vy(i+1) = vy(i) + (1/6)*(k_1vy+2*k_2vy+2*k_3vy+k_4vy)*h;
py(i+1) = py(i) + (1/6)*(k_1py+2*k_2py+2*k_3py+k_4py)*h;
z(i+1) = z(i) + (1/6)*(k_1z+2*k_2z+2*k_3z+k_4z)*h;
vz(i+1) = vz(i) + (1/6)*(k_1vz+2*k_2vz+2*k_3vz+k_4vz)*h;
pz(i+1) = pz(i) + (1/6)*(k_1pz+2*k_2pz+2*k_3pz+k_4pz)*h;
end
1 Comment
James Tursa
on 31 Mar 2023
I would suggest you rewrite this code using vectors for your states instead of individual variables. E.g., define a state vector y as follows:
y(1) = x
y(2) = y
y(3) = z
y(4) = vx
y(5) = vy
y(6) = vz
y(7) = px
y(8) = py
y(9) = pz
Then define one derivative function handle:
f = @(t,y) = [y(4:6); other stuff ]
Your downstream code will be greatly simplified and easier to debug, and this f can be passed directly into ode45( ) for comparison results.
Answers (2)
Cris LaPierre
on 30 Mar 2023
Edited: Cris LaPierre
on 30 Mar 2023
The error is because you define fx with 9 inputs, but call it with 10.
fx = @(x,vx,y,vy,z,vz,px,py,pz) vx;
% vs
k_1x = fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
It looks like you left out time when defining fx. Did you mean to do this (add time to function inputs)?
fx = @(t,x,vx,y,vy,z,vz,px,py,pz) vx;
or this (remove time from calling syntax)
k_1x = fx(x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
It looks like you'll have this same issue with all your functions.
6 Comments
Ban
on 30 Mar 2023
Thanks. Yeah, I missed out the time variable as a typo.
It's too many inputs and too many eqns. I would like to vectorize everything, like this - X=[x,vx,y,vy,z,vx,px,py,pz]; F=@(t,X)[...all the eqns...]; Then vectorize K1,K2,K3,K4.
But the issue is it shows error stating that variables in X are undefined. Please suggest.
Torsten
on 30 Mar 2023
I wonder why you program the numerical method on your own.
MATLAB has built-in efficient and well-tested ODE integrators you can use for your problem.
This way, you can be quite certain that problems you get with the code are due to your equations, not due to the self-programmed Runge-Kutta ODE integrator.
Cris LaPierre
on 30 Mar 2023
Please share the full error message you are seeing (all the red text)
AkB
on 30 Mar 2023
Walter Roberson
on 30 Mar 2023
That line does not exist in your posted code, not even as a comment.
AkB
on 31 Mar 2023
I changed
(1./(1+(1/B).*C22d)).*(-c1.*(K22.*vy));
to
(1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
in the definition of fvy.
And you should check these strange cos(acos(...)) expressions in your differential equations.
ae=10 ; % aspect ratio
ds=1*10^-3 ; % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ; % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2; % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
Tp=2*pi/w1 ; % time period
B=1.1;
theta_ini= 0 ; % azimuthal angle matrix
phi_ini= pi/4 ; % polar angle matrix
vp= 0.1*w1*A ; % Particle intrinsic velocity (m/s)
H=3700; % depth of ocean
Fr=w1*A./sqrt(g*H); % Froude number
z0= 3 ;
c1=18*nu/(B*ds^2) ;
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*Tp; %seconds
dt=0.05 ;
tspan=t0:dt:tf ;
% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;
u0y=0;
u0z=0;
% Initial particle velocity
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;
% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);
x = 0 ;
vx = v_ini_1p ;
y = 0 ;
vy = v_ini_2p ;
z = -z0 ;
vz = v_ini_3p ;
px = initial_p1 ;
py = initial_p2 ;
pz = initial_p3 ;
Y0 = [x;vx;y;vy;z;vz;px;py;pz];
[T,Y] = ode45(@fun,tspan,Y0);
plot(T,Y(:,1))
function dY = fun(t,Y)
% Y=[x; vx; y; vy; z; vz; px; py; pz];
x = Y(1);
vx = Y(2);
y = Y(3);
vy = Y(4);
z = Y(5);
vz = Y(6);
px = Y(7);
py = Y(8);
pz = Y(9);
ae=10 ; % aspect ratio
ds=1*10^-3 ; % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ; % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2; % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ; % time period
B=1.1;
theta_ini= 0 ; % azimuthal angle matrix
phi_ini= pi/4 ; % polar angle matrix
vp= 0.1*w1*A ; % Particle intrinsic velocity (m/s)
H=3700; % depth of ocean
Fr=w1*A./sqrt(g*H); % Froude number
z0= 3 ;
c1=18*nu/(B*ds^2) ;
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
dY(1) = vx;
dY(2) = (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
(1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t))));
dY(3) = vy;
dY(4) = (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
dY(5) = vz;
dY(6) = (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
(1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g);
dY(7) = ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
dY(8) = -ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
dY(9) = ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
dY = dY.';
end
3 Comments
Your code works, too:
ae=10 ; % aspect ratio
ds=1*10^-3 ; % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ; % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2; % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ; % time period
B=1.1;
theta_ini= 0 ; % azimuthal angle matrix
phi_ini= pi/4 ; % polar angle matrix
vp= 0.1*w1*A ; % Particle intrinsic velocity (m/s)
H=3700; % depth of ocean
Fr=w1*A./sqrt(g*H); % Froude number
z0= 3 ;
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*T; %seconds
dt=0.05 ;
t=t0:dt:tf ;
c1=18*nu/(B*ds^2) ;
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
% X=[x; vx; y; vy; z; vz; px; py; pz];
% RK4
fx = @(t,x,vx,y,vy,z,vz,px,py,pz) vx;
fvx = @(t,x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
(1/B)*((C11d.*sin(acos(pz)).*sin(acos(pz))+C22d.*cos(acos(pz)).*cos(acos(pz)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
c1*((K11d.*sin(acos(pz)).*sin(acos(pz))+K22d.*cos(acos(pz)).*cos(acos(pz))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t))));
fy = @(t,x,vx,y,vy,z,vz,px,py,pz) vy;
fvy = @(t,x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
fz = @(t,x,vx,y,vy,z,vz,px,py,pz) vz;
fvz = @(t,x,vx,y,vy,z,vz,px,py,pz) (1./(1+(1/B).*(C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
(1/B)*((C11d.*cos(acos(pz)).*cos(acos(pz))+C22d.*sin(acos(pz)).*sin(acos(pz)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
c1.*((K11d.*cos(acos(pz)).*cos(acos(pz))+K22d.*sin(acos(pz)).*sin(acos(pz))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g);
fpx = @(t,x,vx,y,vy,z,vz,px,py,pz) ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz)-ecc.*px.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpy = @(t,x,vx,y,vy,z,vz,px,py,pz) -ecc.*py.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+...
(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
fpz = @(t,x,vx,y,vy,z,vz,px,py,pz) ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz)-ecc.*pz.*((px.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*px+...
pz.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*px)+(px.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*pz+pz.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*pz));
h=dt;
x=zeros(length(t),1);
vx=zeros(length(t),1);
px=zeros(length(t),1);
y=zeros(length(t),1);
vy=zeros(length(t),1);
py=zeros(length(t),1);
z=zeros(length(t),1);
vz=zeros(length(t),1);
pz=zeros(length(t),1);
% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;
u0y=0;
u0z=0;
% Initial particle velocity
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;
% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);
x(1) = 0 ;
vx(1) = v_ini_1p ;
y(1) = 0 ;
vy(1) = v_ini_2p ;
z(1) = -z0 ;
vz(1) = v_ini_3p ;
px(1) = initial_p1 ;
py(1) = initial_p2 ;
pz(1) = initial_p3 ;
for i=1: length(t)-1
k_1x = fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1vx = fvx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1y = fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1vy = fvy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1z = fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1vz = fvz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
K1y = fy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
K1z = fz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1px = fpx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1py = fpy(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_1pz = fpz(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
k_2x = fx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2y = fy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2z = fz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_2pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_1x, vx(i)+0.5*h*k_1vx, y(i)+0.5*h*k_1y, vy(i)+0.5*h*k_1vy, z(i)+0.5*h*k_1z, vz(i)+0.5*h*k_1vz, px(i)+0.5*h*k_1px, py(i)+0.5*h*k_1py, pz(i)+0.5*h*k_1pz);
k_3x = fx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3y = fy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3z = fz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_3pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_2x, vx(i)+0.5*h*k_2vx, y(i)+0.5*h*k_2y, vy(i)+0.5*h*k_2vy, z(i)+0.5*h*k_2z, vz(i)+0.5*h*k_2vz, px(i)+0.5*h*k_2px, py(i)+0.5*h*k_2py, pz(i)+0.5*h*k_2pz);
k_4x = fx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4vx = fvx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4y = fy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4vy = fvy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4z = fz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4vz = fvz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4px = fpx(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4py = fpy(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
k_4pz = fpz(t(i)+0.5*h, x(i)+0.5*h*k_3x, vx(i)+0.5*h*k_3vx, y(i)+0.5*h*k_3y, vy(i)+0.5*h*k_3vy, z(i)+0.5*h*k_3z, vz(i)+0.5*h*k_3vz, px(i)+0.5*h*k_3px, py(i)+0.5*h*k_3py, pz(i)+0.5*h*k_3pz);
x(i+1) = x(i) + (1/6)*(k_1x+2*k_2x+2*k_3x+k_4x)*h;
vx(i+1) = vx(i) + (1/6)*(k_1vx+2*k_2vx+2*k_3vx+k_4vx)*h;
px(i+1) = px(i) + (1/6)*(k_1px+2*k_2px+2*k_3px+k_4px)*h;
y(i+1) = y(i) + (1/6)*(k_1y+2*k_2y+2*k_3y+k_4y)*h;
vy(i+1) = vy(i) + (1/6)*(k_1vy+2*k_2vy+2*k_3vy+k_4vy)*h;
py(i+1) = py(i) + (1/6)*(k_1py+2*k_2py+2*k_3py+k_4py)*h;
z(i+1) = z(i) + (1/6)*(k_1z+2*k_2z+2*k_3z+k_4z)*h;
vz(i+1) = vz(i) + (1/6)*(k_1vz+2*k_2vz+2*k_3vz+k_4vz)*h;
pz(i+1) = pz(i) + (1/6)*(k_1pz+2*k_2pz+2*k_3pz+k_4pz)*h;
end
plot(t,x)
Just a query: when I used the same eqns with another different code, it gives different results. Please check.
tic
clear;
close all;
clc;
ae=10 ; % aspect ratio
ds=1*10^-3 ; % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ; % eccentricity
g=9.81;
nu=0.89*10^-6 ;
A=0.2; % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ; % time period
B=1.1;
theta_ini= 0 ; % azimuthal angle matrix
phi_ini= pi/4 ; % polar angle matrix
vp= 0.1*w1*A ; % Particle intrinsic velocity (m/s)
theta_ini_matrix =[0 pi/4 pi/2] ; % azimuthal angle matrix
phi_ini_matrix =[0 pi/4 pi/2] ; % polar angle matrix
vp_matrix=[10^-1 10^-2 10^-3].*w1*A ; % Particle intrinsic velocity (m/s)
H=3700; % depth of ocean
Fr=w1*A./sqrt(g*H); % Froude number
% Initial depth of particles
z0= 3 ;
% time scale of photoaxis
tau_s_matrix =[0.1 1 10] ;
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20*T; %seconds
interval = [t0 tf]; %time interval
n_points=200;
tValues = linspace(interval(1),interval(2),tf*n_points);
syms p1(t) p2(t) p3(t) x(t) y(t) z(t)
% mass tensor in body axes
alpha0=(1./((ae.^2)-1).^(3/2)).*(2*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
beta0=(0.5./((ae.^2)-1).^(3/2)).*(2.*(ae.^2).*sqrt((ae.^2)-1)+ae.*log((ae-sqrt(ae.^2-1))./(ae+sqrt(ae.^2-1)))) ;
gamma0=beta0 ;
C11d=alpha0./(2-alpha0) ;
C22d=beta0./(2-beta0) ;
C33d=gamma0./(2-gamma0) ;
% resistance tensor in body axes
K11d=((4/3).*(ae.^(-1/3)).*(1-ae.^2))./(ae-((2*(ae.^2)-1).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K22d=((8/3)*(ae.^(-1/3)).*((ae.^2)-1))./(ae+((2*(ae.^2)-3).*log(ae+sqrt(ae.^2-1))./sqrt(ae.^2-1))) ;
K33d=K22d ;
K33_ini=K11d*cos(phi_ini)*cos(phi_ini)+K22d*sin(phi_ini)*sin(phi_ini) ; % K33 @t=0
% Defining the governing equations
c1=18*nu/(B*ds^2) ;
vx=diff(x,t);
vy=diff(y,t);
vz=diff(z,t);
% Maxey–Riley Equation (passive)
eq1 = diff(x,t,t) == (1./(1+(1/B).*(C11d.*sin(acos(p3)).*sin(acos(p3))+C22d.*cos(acos(p3)).*cos(acos(p3))))).*((1/B).*((w1^2)*A*exp(k*z)*sin(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(-k*w1*A*exp(k*z)*sin(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t)))+...
(1/B)*((C11d.*sin(acos(p3)).*sin(acos(p3))+C22d.*cos(acos(p3)).*cos(acos(p3)))*(w1^2)*A*exp(k*z)*sin(k*x-w1*t))-...
c1*((K11d.*sin(acos(p3)).*sin(acos(p3))+K22d.*cos(acos(p3)).*cos(acos(p3))).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t))));
eq2 = diff(y,t,t) == (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));
eq3 = diff(z,t,t) == (1./(1+(1/B).*(C11d.*cos(acos(p3)).*cos(acos(p3))+C22d.*sin(acos(p3)).*sin(acos(p3))))).*((1/B).*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)+...
w1*A*exp(k*z)*cos(k*x-w1*t).*(k*w1*A*exp(k*z)*cos(k*x-w1*t))+w1*A*exp(k*z)*sin(k*x-w1*t).*(k*w1*A*exp(k*z)*sin(k*x-w1*t)))+...
(1/B)*((C11d.*cos(acos(p3)).*cos(acos(p3))+C22d.*sin(acos(p3)).*sin(acos(p3)))*(-(w1^2)*A*exp(k*z)*cos(k*x-w1*t)))-...
c1.*((K11d.*cos(acos(p3)).*cos(acos(p3))+K22d.*sin(acos(p3)).*sin(acos(p3))).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g);
% Jeffery’s Equation
eq4 = diff(p1,t) == ecc.*((-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3)-ecc.*p1.*((p1.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+...
p3.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1)+(p1.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3+p3.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3));
eq5 = diff(p2,t) == -ecc.*p2.*((p1.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+p3.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1)+...
(p1.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3+p3.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3));
eq6 = diff(p3,t) == ecc.*((k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1+(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3)-ecc.*p3.*((p1.*(-k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p1+...
p3.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p1)+(p1.*(k*w1*A*exp(k*z)*cos(k*x-w1*t)).*p3+p3.*(k*w1*A*exp(k*z)*sin(k*x-w1*t)).*p3));
vars = [x(t); y(t); z(t); p1(t); p2(t); p3(t)];
V_p = odeToVectorField([eq2 eq1 eq3 eq4 eq5 eq6])
M_p = matlabFunction(V_p,'vars', {'t','Y'});
% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;
u0y=0;
u0z=0;
% Initial particle velocity
v_ini_1p=u0x ;
v_ini_2p=u0y ;
v_ini_3p=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;
% Initial particle orientation
initial_p1 = sin(phi_ini).*cos(theta_ini);
initial_p2 = sin(phi_ini).*sin(theta_ini);
initial_p3 = cos(phi_ini);
%# Solving for passive particle
y0_p=[0 v_ini_1p 0 v_ini_2p 0 v_ini_3p initial_p1 initial_p2 initial_p3]; %Initial Conditions for passive particle
ySol_p = ode45(M_p,interval,y0_p); %solution for passive particle
x_passive = deval(ySol_p,tValues,1);
y_passive = deval(ySol_p,tValues,3);
z_passive = -z0+deval(ySol_p,tValues,5);
v1_passive = deval(ySol_p,tValues,2);
v2_passive = deval(ySol_p,tValues,4);
v3_passive = deval(ySol_p,tValues,6);
p1_passive = deval(ySol_p,tValues,7);
p2_passive = deval(ySol_p,tValues,8);
p3_passive = deval(ySol_p,tValues,9);
phi_passive=acos(p3_passive) ;
theta_passive=asin(real(p2_passive./sin(phi_passive)));
%# plotting of results
plot(tValues,x_passive,'LineWidth',2)
xlabel('$t$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$x$', 'FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);
toc
Categories
Find more on Programming in Help Center and File Exchange
on 31 Mar 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
