Runge Kutta 4th Order Method

조회 수: 22(최근 30일)
Ryan Bowman
Ryan Bowman 2020년 3월 4일
답변: Meysam Mahooti 2021년 5월 5일
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.)
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)
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
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일
Also, in your function, in the statenent
Imat_inv = diag(1./I);
the inner operation
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:
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일

Community Treasure Hunt

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

Start Hunting!

Translated by