Runge-kutta - free-fall motion

7 次查看（过去 30 天）
Noel Lou
Noel Lou 2016-3-5
编辑： John D'Errico 2016-3-5
Hi there, i am doing a trajectory simulation for a free-fall lifeboat, however i tried solving the following motion equations to produce a trajectory as shown in the image attached
X' = Vx
Z' = Vz
Vx'=[fn*(sin θ - uCos θ )]/M
Vz'={[fn(cos θ + uSin θ )]/M} - g
however my runge kutta code produce something different, could someone help me see where's my mistake?
thank you
clc; % Clears the screen
clear all;
thete=30;
g= 9.81; %Constant
h=0.001; % step size
tfinal = 3;
N=ceil(tfinal/h); %ceil is round up
t(1) = 0;% initial condition
Vx(1)=0;%initial accleration
X(1)=0;
Vz(1)=0;
Z(1)=0;%initial velocity
F_X = @(t,X,Vx) Vx;
F_Z = @(t,Z,Vz) Vz;
F_Vx = @(t,X,Vx)(0.866*(sin(thete)-0.5774*(cos(thete))));
F_Vz = @(t,Z,Vz)(0.866*(cos(thete)+0.5774*(sin(thete)))-g);
for i=1:N-1; % calculation loop
%update time
t(i+1)=t(i)+h;
%update motions main equation
k_1 = F_X (t(i) ,X(i) ,Vx(i));
L_1 = F_Vx(t(i) ,X(i) ,Vx(i));
k_2 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1);
L_2 = F_Vx(t(i)+0.5*h,X(i)+0.5*h*k_1,Vx(i)+0.5*h*L_1);
k_3 = F_X (t(i)+0.5*h,X(i)+0.5*h*k_2,Vx(i)+0.5*h*L_2);
L_3 = F_Vx(t(i)+0.5*h,X(i)+0.5*h*k_2,Vx(i)+0.5*h*L_2);
k_4 = F_X (t(i)+h ,X(i)+h *k_3,Vx(i)+ h*L_3);
L_4 = F_Vx(t(i)+h ,X(i)+h *k_3,Vx(i)+ h*L_3);
k_11 = F_Z(t(i) ,Z(i) ,Vz(i));
L_11 = F_Vz(t(i) ,Z(i) ,Vz(i));
k_22 = F_Z(t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11);
L_22 = F_Vz(t(i)+0.5*h,Z(i)+0.5*h*k_11,Vz(i)+0.5*h*L_11);
k_33 = F_Z(t(i)+0.5*h,Z(i)+0.5*h*k_22,Vz(i)+0.5*h*L_22);
L_33 = F_Vz(t(i)+0.5*h,Z(i)+0.5*h*k_22,Vz(i)+0.5*h*L_22);
k_44 = F_Z(t(i)+ h,Z(i)+ h*k_33,Vz(i)+ h*L_33);
L_44 = F_Vz(t(i)+ h,Z(i)+ h*k_33,Vz(i)+ h*L_33);
X(i+1) = X(i) + (h/6)*(k_1+2*k_2+2*k_3+k_4);
Vx(i+1) = X(i) + (h/6)*(L_1+2*L_2+2*L_3+L_4);
Z(i+1) = Z(i) + (h/6)*(k_11+2*k_22+2*k_33+k_44);
Vz(i+1) = Z(i) + (h/6)*(L_11+2*L_22+2*L_33+L_44);
end
plot(t,X,'--', 'Linewidth', 1.5, 'color', 'red')
hold on
plot(t,Vx,'--', 'Linewidth', 1.5, 'color', 'black')
hold on
plot(t,Z,'--', 'Linewidth', 1.5, 'color', 'blue')
hold on
plot(t,Vz,'--', 'Linewidth', 1.5, 'color', 'yellow')
xlabel('time')
ylabel('height')
legend('X','Vx','Z','Vz')
  3 个评论
显示 1更早的评论
Noel Lou
Noel Lou 2016-3-5
hi there,
thank so much for your reply, however i tried ode45 and its not a stiff system, it produce the same graph as my above RK4, however i notice that at the right handside of my function all i is got a constant , and does not have any x,z,vx,vz,t .
t(1) = 0;% initial condition
Vx(1)=0;%initial accleration
X(1)=0;
Vz(1)=0;
Z(1)=0;%initial velocity
F_X = @(t,X,Vx) Vx*t;
F_Z = @(t,Z,Vz) Vz*t;
F_Vx = @(t,X,Vx)(0.866*(sin(thete)-0.5774*(cos(thete))))*t;
F_Vz = @(t,Z,Vz)(0.866*(cos(thete)+0.5774*(sin(thete)))-9.81)*t;
could it be due to this which cause my graph to look not like the picture above, im trying to simulate the sliding phase, so my graph should look like the slope at the top left hand side of the picture up till 5m in x-position.
thank you
John D'Errico
John D'Errico 2016-3-5
编辑：John D'Errico 2016-3-5
Perhaps the first question one should always ask about an ODE problem is if the system is stiff. However, my point was that ODE45 will also fail on stiff problems, while a solver designed for that purpose, like ode15s or ode23s, will succeed. Regardless, I will guess the issue is NOT a stiffness one here.
So in order for me to understand your code, I need to understand how you are modeling the system. You have a falling "object", that apparently starts out above the water, and eventually hits the water and goes under a bit, then pops up?
From the expected plot, it appears that there is vertical motion,as well as horizontal motion.
t is clearly time
X is apparently position, so we start out at 0. Can I assume that X(t) is the position of the center of the object?
Z is also position, apparently height above the water?
thete is apparently an angle, so it appears the object is tilted at an angle of thete initially?
You use variables Vx and Vz, apparently as velocities. But there appears to be a horizontal component of velocity at the start in that plot. That belies the statement that Vx is zero at the start. As well, I don't see anything in the model that would cause the object to move horizontally, UNTIL it would hit the water.
h appears to be a step size, so a time step increment. (A more mnemonic variable name might have been dt. Using easy to read and follow variable names is very helpful in trying to read code.)
tfinal is 3, so I assume you will simulate the process for 3 seconds.
There are some other variables in the equation as provided on top, that are not in your code. For example, M is probably the mass of the object, and G is clearly the gravitational constant. u and fn seem to be missing completely, though you seem to have put in some numbers there.
Sorry for asking these questions, but when I read your code, and try to compare it to the comments inline and what I've seen so far, there are inconsistencies.
For example, logically, we are told that the initial position is 27m ABOVE the water. But Z(1) is set to 0. More confusing is this line:
Z(1)=0;%initial velocity
which seems to indicate that Z is a velocity. However, it would appear that Vz should be a Z direction velocity. I would also have expected that Z(1) would be 27, so the height above the water.
My guess about your problem is that you are incorrectly transcribing the mathematics into the equations of motion, but until I understand the system you are modeling, it is difficult to help you.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by