Variable input ODE solution
1 次查看(过去 30 天)
显示 更早的评论
The following code is for solving a system of ODEs for a fixed input B (corresponding to jB=1).
I would like to modify it to a case where, as time (t) starts, the system of ODEs are solved for B (corresponding to jB=1). However, when vx becomes vx <0, then the system of ODEs should be solved for B (corresponding to jB=2). Again, when vx becomes vx <0 for jB=2 case, then the system of ODEs should be again solved for B (corresponding to jB=1). And continue the process through repeatations. It is to be noted that, every time jB values are switched, x, y, z, px, py, pz should be continued from the earlier case. This is because, it is a particle that moves with jB=1 untill vx <0, then jB=2 and the particles sinks from the position at which vx <0.
Please suggest how to proceed.
My code is:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;close all; clc;
nu=0.89*10^-6 ;
g=9.81;
A=0.2; % wave amplitude
k=0.33/A; % k*A <= 0.33 for Linear wave theory to hold correct
L=2*pi/k ;
w1=sqrt(g*k); % angular frequency
T=2*pi/w1 ; % time period
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =10*T; %seconds
tspan=[t0 tf] ;
h=0.01;
ecc_matrix=[0.2 0.5 0.9]; % eccentricity
theta_ini_matrix =[0 pi/4 pi/2] ; % azimuthal angle matrix
phi_ini_matrix =[0 pi/4 pi/2] ; % polar angle matrix
B_matrix=[1.0001 5]; % B=particle density/fluid density matrix
ds_matrix=[1*10^-3 1*10^-4]; % equivalent dia matrix
x0=0; y0=0; z0=-0.1;
for jds=1
for jtheta=1
for jphi=1
ds=ds_matrix(jds) ; % equivalent dia
theta_ini = theta_ini_matrix(jtheta) ; % azimuthal angle
phi_ini = phi_ini_matrix(jphi) ; % polar angle
% Initial field velocity @t=0, x=0, z=0
u0x=w1*A*exp(k*z0)*cos(k*x0);
u0y=0;
u0z=w1*A*exp(k*z0)*sin(k*x0);
% Initial particle orientation
p0x = sin(phi_ini).*cos(theta_ini);
p0y = sin(phi_ini).*sin(theta_ini);
p0z = cos(phi_ini);
for jecc=2
ecc=ecc_matrix(jecc);
ae=sqrt((1+ecc)./(1-ecc));
% resistance tensor in body axes (prolate)
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
% 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) ;
t = tspan(1):h:tspan(2);
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);
u0x=w1*A*exp(k*z0)*cos(k*x0);
u0y=0;
u0z=w1*A*exp(k*z0)*sin(k*x0);
for i=1: length(t)-1
for jB=1
B=B_matrix(jB) ; % B=particle density/fluid density
c1=18*nu/(B*ds^2) ;
Fx = @(t,x,vx,y,vy,z,vz,px,py,pz) [vx; (1./(1+(1/B).*(C11d.*(px.*px+py.*py)+C22d.*pz.*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.*(px.*px+py.*py)+C22d.*pz.*pz)*((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)))-c1*((K11d.*(px.*px+py.*py)+K22d.*pz.*pz).*(vx-w1*A*exp(k*z)*cos(k*x-w1*t))));vy;...
(1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));vz;(1./(1+(1/B).*(C11d.*pz.*pz+C22d.*(px.*px+py.*py)))).*((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.*pz.*pz+...
C22d.*(px.*px+py.*py))*((-(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))))-c1.*((K11d.*pz.*pz+K22d.*(px.*px+py.*py)).*(vz-w1*A*exp(k*z)*sin(k*x-w1*t)))-(1-1/B).*g);...
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));...
-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)); 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))];
v0x=u0x ;
v0y=u0y ;
v0z=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z ;
x(1) = x0 ;
vx(1) = v0x ;
y(1) = y0 ;
vy(1) = v0y ;
z(1) = z0 ;
vz(1) = v0z ;
px(1) = p0x ;
py(1) = p0y ;
pz(1) = p0z ;
K1 = Fx(t(i),x(i),vx(i),y(i),vy(i),z(i),vz(i),px(i),py(i),pz(i));
K2 = Fx(t(i)+0.5*h, x(i)+0.5*h*K1(1), vx(i)+0.5*h*K1(2), y(i)+0.5*h*K1(3), vy(i)+...
0.5*h*K1(4), z(i)+0.5*h*K1(5), vz(i)+0.5*h*K1(6), px(i)+0.5*h*K1(7), py(i)+0.5*h*K1(8), pz(i)+0.5*h*K1(9));
K3 = Fx(t(i)+0.5*h, x(i)+0.5*h*K2(1), vx(i)+0.5*h*K2(2), y(i)+0.5*h*K2(3), vy(i)+...
0.5*h*K2(4), z(i)+0.5*h*K2(5), vz(i)+0.5*h*K2(6), px(i)+0.5*h*K2(7), py(i)+0.5*h*K2(8), pz(i)+0.5*h*K2(9));
K4 = Fx(t(i)+h, x(i)+h*K3(1), vx(i)+h*K3(2), y(i)+h*K3(3), vy(i)+h*K3(4), z(i)+...
h*K3(5), vz(i)+h*K3(6), px(i)+h*K3(7), py(i)+h*K3(8), pz(i)+h*K3(9));
x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+K4(1))*h;
vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;
y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+K4(3))*h;
vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+K4(5))*h;
vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
px(i+1) = px(i) + (1/6)*(K1(7)+2*K2(7)+2*K3(7)+K4(7))*h;
py(i+1) = py(i) + (1/6)*(K1(8)+2*K2(8)+2*K3(8)+K4(8))*h;
pz(i+1) = pz(i) + (1/6)*(K1(9)+2*K2(9)+2*K3(9)+K4(9))*h;
end
end
x_passive=x;
y_passive=y;
z_passive=z ;
px_passive=px;
py_passive=py;
pz_passive=pz;
vx_passive=vx;
vy_passive=vy;
vz_passive=vz;
% % plotting of results
plot(k*x_passive,k*z_passive,'-r','LineWidth',2)
xlabel('$kx$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
ylabel('$kz$','FontSize',20,'FontWeight','bold', 'Interpreter','latex');
set(gca,'FontSize',15);
hold on
end
end
end
end
3 个评论
Torsten
2023-7-25
You won't be successful with your own Runge-Kutta code to accomplish this.
Use ode45 and its event locating facility.
回答(1 个)
Mrutyunjaya Hiremath
2023-7-25
Here It is ..
clear; close all; clc;
% (Please insert the rest of your code
for jds = 1:length(ds_matrix)
for jtheta = 1:length(theta_ini_matrix)
for jphi = 1:length(phi_ini_matrix)
ds = ds_matrix(jds); % equivalent dia
theta_ini = theta_ini_matrix(jtheta); % azimuthal angle
phi_ini = phi_ini_matrix(jphi); % polar angle
% ... Rest of your code to initialize other parameters ...
for jecc = 2
ecc = ecc_matrix(jecc);
ae = sqrt((1 + ecc) / (1 - ecc));
% ... Rest of your code to initialize resistance tensor and mass tensor ...
t = tspan(1):h:tspan(2);
x = zeros(length(t), 1);
vx = zeros(length(t), 1);
% ... Initialize other variables (y, vy, z, vz, px, py, pz) ...
% Set initial jB value to 1
jB = 1;
% Loop through each time step
for i = 1:length(t) - 1
% Check if vx is negative to switch jB
if vx(i) < 0
jB = 3 - jB; % Toggle between 1 and 2 (if jB=1, then jB=2; if jB=2, then jB=1)
end
%... Continue solving the ODEs based on the current jB value ...
% (Please insert the rest of your code for solving ODEs here)
end
%... Continue with the rest of your code ...
end
end
end
end
% ... Rest of your code ...
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!