Runge Kutta 4th Order Method

5 次查看（过去 30 天）
Ryan Bowman
Ryan Bowman 2020-3-4
回答： gn 2023-12-13
I have developed a 4th order runge kutta method that helps me find angular velocity of a rigid body. It should be displayed as a 3x1 matrix to include x,y and z velocites. Unfornately, my k2,k3, and k4 are 3x3 matrices and my k1 is a 3x1 matrix. I believe all matrices should be a 3x1. Here is my code: (I iuncluded my function which is a seperate file.)
clc
clear all
clear variables
format short e
%% Problem 1
I = [5.3 0 0;0 7.4 0;0 0 10.5];
W_0 = [0.1 0.2 0.3];
T = [0 0 0]'; %No torque is given in this case
dt = 0.025; %seconds
t = 0:dt:10; %seconds
[Wdot] = getWdot(W_0,T,I);
W(1,:) = W_0;
for ii = 1:length(t)-1
k1 = dt*getWdot(W(ii,:),T,I);
k2 = dt*getWdot(W(ii,:)+0.5*k1, t+0.5*dt, I);
k3 = dt*getWdot(W(ii,:)+0.5*k2, t+0.5*dt, I);
k4 = dt*getWdot(W(ii,:)+k3, t+dt, I);
W(ii+1,:) = W(ii,:) + (1/6)*(k1+2*k2+2*k3+k4)
end
function [Wdot] = getWdot(W_0,T,I)
Imat = diag(I);
Imat_inv = diag(1./I);
Wdot = Imat_inv.*cross(W_0',Imat.*W_0')
Any suggestions to fix this issue would be much appreciated.

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

回答（4 个）

Jim Riggs
Jim Riggs 2020-3-5
编辑：Jim Riggs 2020-3-5
I think the problem is the T argument.
in k1 you have a capital T, which is a 3-vector
in k2, k3, and k4 it is a lower case t, which is 0:dt:10.
Something is wrong there.
k1 = dt*getWdot(W(ii,:), T, I);
k2 = dt*getWdot(W(ii,:)+0.5*k1, t+0.5*dt, I);
k3 = dt*getWdot(W(ii,:)+0.5*k2, t+0.5*dt, I);
k4 = dt*getWdot(W(ii,:)+k3, t+dt, I);
They are probably both wrong. I think what you want is t(ii) (although I don't see t used in the function)
  1 个评论
Jim Riggs
Jim Riggs 2020-3-5
编辑：Jim Riggs 2020-3-5
Also, in your function, in the statenent
Imat_inv = diag(1./I);
the inner operation
(1./I)
produces a 3x3 matrix wth divide by zero in the off diagonal terms.
You would be better off using
Imat_inv = 1./Imat;
since Imat is already a 3-vector. There are no divide by zeros in this operatin and you end up with the same result.
Otherwise, you should use
Imat_inv = diag(inv(I));
This produces a propper matrix inverse, from which you can select the diagonal.

请先登录，再进行评论。

James Tursa
James Tursa 2020-3-5
编辑：James Tursa 2020-3-5
When solving the EOM for w_dot you should get a minus sign:
https://en.wikipedia.org/wiki/Euler%27s_equations_%28rigid_body_dynamics%29
I.e., this
Wdot = Imat_inv.*cross(W_0',Imat.*W_0')
should be this instead:
Wdot = -Imat_inv.*cross(W_0',Imat.*W_0')
Also, I guess you are trying to be clever in calculating the inertia matrix inverse and multiplies with all of that diag stuff. Why bother? Just use backslash:
function [Wdot] = getWdot(W_0,T,I)
Wdot = -I \ cross(W_0',I*W_0')
And, you are inconsistent in your third argument T above. For k1 you pass in T, but in k2, k3, k4 you pass in an entire time vector. You don't use this argument in the derivative so it didn't cause any problems, but your calling statements are incorrect ... you probably mean to use t(ii) instead of t for that argument.
Finally, although not necessary, I would have opted to have all your vectors be column vectors so that you didn't need to put in all of these transposes for the calculations. But that is just a nit.
Meysam Mahooti
Meysam Mahooti 2021-5-5
https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle
gn
gn 2023-12-13
x = x0
y = y0
x_values = [x0]
y_values = [y0]
while x < xn:
k1 = h * dy_dx(x, y)
k2 = h * dy_dx(x + 0.5 * h, y + 0.5 * k1)
k3 = h * dy_dx(x + 0.5 * h, y + 0.5 * k2)
k4 = h * dy_dx(x + h, y + k3)
y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0
x = x + h
x_values.append(x)
y_values.append(y)
return x_values, y_values
# Example usage:
def example_derivative(x, y):
return x + y
x0 = 0
xn = 1
y0 = 1
h = 0.1
x_values, y_values = runge_kutta_4th_order(example_derivative, x0, y0, xn, h)

类别

在 Help Center 和 File Exchange 中查找有关 Runge Kutta Methods 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by