Problem with time steps in Runge-Kutte method
2 次查看(过去 30 天)
显示 更早的评论
I am solving 6 ODEs with Runge-Kutte method. But I get different results with different time steps. Please suggest how to deal with this. My code is following:
clear;
close all;
clc;
nu=0.89*10^-6 ;
g=9.81;
ae=10 ; % aspect ratio
ds=1*10^-3 ; % equivalent dia
ecc=((ae.^2)-1)./((ae.^2)+1) ; % eccentricity
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
zh= 3 ;
vp= 0.01*w1*A ; % Particle intrinsic velocity (m/s)
% Time span for evolution of dynamics
t0 = 0; %seconds
tf =20; %seconds
tspan=[t0 tf] ;
tau_s=5 ;
h_mat=[0.1 0.05 0.01 0.005];
for jh=1:4
h=h_mat(jh);
% 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
x0=0;
y0=0;
z0=0;
% Initial field velocity @t=0, x=0, z=0
u0x=w1*A;
u0y=0;
u0z=0;
% Initial particle orientation
p0x = sin(phi_ini).*cos(theta_ini);
p0y = sin(phi_ini).*sin(theta_ini);
p0z = cos(phi_ini);
% Initial particle velocity
v0x=u0x+vp*p0x ;
v0y=u0y+vp*p0y ;
v0z=-((B-1)*g*ds^2)/(18*nu*K33_ini)+u0z+vp*p0z ;
nu=0.89*10^-6 ;
c1=18*nu/(B*ds^2) ;
g=9.81;
% 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);
dx=zeros(length(t),1);
dvx=zeros(length(t),1);
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 ;
Fx = @(t,x,vx,y,vy,z,vz,px,py,pz) [vx; (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)))); vy; (1./(1+(1/B).*C22d)).*(-c1.*(K22d.*vy));vz; (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); -px.*pz/tau_s+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));-py.*pz/tau_s-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));(1-pz.*pz)/tau_s+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))];
for i=1: length(t)-1
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));
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;
vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;
vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
vx(i+1) = vx(i+1) + vp*px(i+1) ;
vy(i+1) = vy(i+1) + vp*py(i+1) ;
vz(i+1) = vz(i+1) + vp*pz(i+1) ;
x(i+1) = x(i)+vx(i+1)*h;
y(i+1) = y(i)+vy(i+1)*h;
z(i+1) = z(i)+vz(i+1)*h;
end
% % plotting of results
plot(k*x, k*z, '-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
0 个评论
回答(2 个)
James Tursa
2023-4-3
编辑:James Tursa
2023-4-3
"... I get different results with different time steps ..."
In general, that is the expected result. I suppose the real question is are any of your results valid? To answer that, I suggest you rewrite your derivative equations in a vector format so that it can be fed to one of the MATLAB solvers such as ode45( ) for comparison. E.g., create this function handle to use with ode45( ):
F = @(t,y) Fx(t,y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8),y(9))
4 个评论
Ban
2023-4-6
My eqns are pretty messy. So in terms of marix it looks very confusing. I am trying to write somrthing like this:
fx = @(t,x,vx,y,vy,z,vz) C11d.*sin(acos(pz)) ;
fy = @(t,x,vx,y,vy,z,vz) C22d.*(K22d.*vy)) ;
fz = @(t,x,vx,y,vy,z,vz) C11d.*cos(acos(pz)) ;
Fx = @(t,x,vx,y,vy,z,vz) [vx; fx; vy; fy; vz; fz; 0; 0; 0];
It gives this error:
The following error occurred converting from double to function_handle:
Too many output arguments.
Torsten
2023-4-6
编辑:Torsten
2023-4-6
Write a function F in which you compute the time derivatives of your solution variables and call this function by the Runge-Kutta solver. Just look at the code I linked to. Your equations don't need to have a matrix structure for this:
function dZ = F(t,Z)
x = Z(1);
vx = Z(2);
y = Z(3);
vy = Z(4);
z = Z(5);
vz = Z(6);
px = Z(7);
py = Z(8);
pz = Z(9);
% Compute time derivatives dx, dvx, dy, dvy, dz, dvz, dpx, dpy dpz without using function handles
...
% Return the vector of time derivatives [dx dvx dy dvy dz dvz dpx dpy dpz] to runge-kutta code
dZ = [dx dvx dy dvy dz dvz dpx dpy dpz].'
end
But I still don't understand why you don't use ODE45. Do you think your Runge-Kutta solver will outperform MATLAB software ?
Torsten
2023-4-3
编辑:Torsten
2023-4-3
vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;
vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
vx(i+1) = vx(i+1) + vp*px(i+1) ;
vy(i+1) = vy(i+1) + vp*py(i+1) ;
vz(i+1) = vz(i+1) + vp*pz(i+1) ;
x(i+1) = x(i)+vx(i+1)*h;
y(i+1) = y(i)+vy(i+1)*h;
z(i+1) = z(i)+vz(i+1)*h;
This is not RK4 for vx, vy, vz, x, y and z.
If you use the usual update, you get the same results independent of h:
vx(i+1) = vx(i) + (1/6)*(K1(2)+2*K2(2)+2*K3(2)+K4(2))*h;
vy(i+1) = vy(i) + (1/6)*(K1(4)+2*K2(4)+2*K3(4)+K4(4))*h;
vz(i+1) = vz(i) + (1/6)*(K1(6)+2*K2(6)+2*K3(6)+K4(6))*h;
x(i+1) = x(i) + (1/6)*(K1(1)+2*K2(1)+2*K3(1)+K4(1))*h;
y(i+1) = y(i) + (1/6)*(K1(3)+2*K2(3)+2*K3(3)+K4(3))*h;
z(i+1) = z(i) + (1/6)*(K1(5)+2*K2(5)+2*K3(5)+K4(5))*h;
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!