Problem solving simultaneous ODEs using 4th order Runge-Kutta method (Too many input arguments.)
古いコメントを表示
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 件のコメント
James Tursa
2023 年 3 月 31 日
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.
回答 (2 件)
Cris LaPierre
2023 年 3 月 30 日
編集済み: Cris LaPierre
2023 年 3 月 30 日
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 件のコメント
Ban
2023 年 3 月 30 日
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
2023 年 3 月 30 日
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
2023 年 3 月 30 日
Please share the full error message you are seeing (all the red text)
AkB
2023 年 3 月 30 日
Walter Roberson
2023 年 3 月 30 日
That line does not exist in your posted code, not even as a comment.
AkB
2023 年 3 月 31 日
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 件のコメント
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
カテゴリ
ヘルプ センター および File Exchange で Programming についてさらに検索
2023 年 3 月 31 日
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
