Variable input ODE solution

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

[Star Strider]

Use the approach described in ODE Event Location to switch between the conditions.

[Ban]

I tried. If you show little more explicitly.

[Torsten]

You won't be successful with your own Runge-Kutta code to accomplish this.

Use ode45 and its event locating facility.

Answer 1

[Mrutyunjaya Hiremath]

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 ...