# Runge Kutta 4th Order Method

38 views (last 30 days)
Ryan Bowman on 4 Mar 2020
Answered: Meysam Mahooti on 5 May 2021
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.

Jim Riggs on 5 Mar 2020
Edited: Jim Riggs on 5 Mar 2020
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 on 5 Mar 2020
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 on 5 Mar 2020
Edited: James Tursa on 5 Mar 2020
When solving the EOM for w_dot you should get a minus sign:
I.e., this
Wdot = Imat_inv.*cross(W_0',Imat.*W_0')
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 on 5 May 2021
https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

R2019a

### Community Treasure Hunt

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

Start Hunting!

Translated by