Runge Kutta 4th Order Method

Question

Ryan Bowman 2020-3-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509080-runge-kutta-4th-order-method

回答： 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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Jim Riggs 2020-3-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509080-runge-kutta-4th-order-method#answer_418588

编辑：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 个评论
显示 -1更早的评论隐藏 -1更早的评论

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.

请先登录，再进行评论。

Answer 2

James Tursa 2020-3-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509080-runge-kutta-4th-order-method#answer_418593

编辑：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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 3

Meysam Mahooti 2021-5-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509080-runge-kutta-4th-order-method#answer_692620

https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 4

gn 2023-12-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509080-runge-kutta-4th-order-method#answer_1370814

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)