Runge Kutta 4th Order Method

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.

Answer 1

[Jim Riggs]

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)

[Jim Riggs]

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]

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.

Answer 3

[Meysam Mahooti]

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

Answer 4

[gn]

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)