2

I am trying to solve a system with seven differential equations. And I am having a hard time grapsing the ode45 solver.
These are the equations:

ω2_dot = -0.75 ω1 ω3

ω1_dot = 0.75 ω2 ω3 + 0.2

ω3_dot = 0

q1_dot = 1/2(ω1q4 + ω2q3 - ω3q2)

q2_dot = 1/2(ω2q4 + ω3q1 - ω1q3)

q3_dot = 1/2(ω3q4 + ω1q2 - ω2q2)

q4_dot = -1/2(ω1q1 + ω2q2 + ω3q3)

The initial values are listed in the inital_val in the same order. Here is what I have so far:

inital_val = [1 -1 2 0 0 0 1];

timespan = [0:0.05:3.95];

[result t] = ode45(@soe,timespan,inital_val);

function [results,t]=soe(inital_val, timespan)
omega1_dot = -0.75*omega2*omega3;
omega2_dot = 0.75*omega2*omega3+0.2;
omega3_dot = 0;
q1_dot = (1/2)*(q4*omega1+omega2*q3-omega3*q2);
q2_dot = (1/2)*(q4*omega2+omega3*q1-omega1*q3);
q3_dot = (1/2)*(q4*omega3+omega1*q2-omega2*q2);
q4_dot = (1/2)*(q1*omega1+q2*omega2+q3*omega3);
results = [omega1 omega2 omega3 q1 q2 q3 q4];

end

However, it gives me this error message, which makes sense but I then don't know how to fix it:

Unrecognized function or variable 'omega2'.

Error in ode45_part>soe (line 10)
omega1_dot = -0.75*omega2*omega3;

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in ode45_part (line 7)
[result t] = ode45(@soe,timespan,inital_val);

Any help would be much appreciated

SHiSA
  • 45
  • 5
  • 1
    Well, you can start by fixing some of the problems with your code. What are the omega and q values in the function? Your function takes `inital_val` as an input, but unless your omegas are defined else where globally, your script does not know how to deal with them. You can. e.g., use something like `[omega1, omega2, omega3, q1, q2, q3, q4] = feval(@(x) x{:}, num2cell(inital_val));` so that it is of the form `omega1 = inital_val(1); omega2 = inital_val(2) ...`, or use something to your liking. – ofey Mar 11 '22 at 19:50
  • You need to take care that the implemented equations correspond to the model equations. The second equation should have w1, not w2. Also check the order of the equations, in the formula section the first two are switched. The last equation has a minus sign as expected from `0.5·q·w` where the squares of the imaginary units are -1 in the real component of the result. – Lutz Lehmann Mar 12 '22 at 07:08

2 Answers2

2

I am not sure what your inital_val is supposed to represent, but here you can at least run the following snippet and fix/change accordingly.

clear; clc;

y0 = [1 -1 2 0 0 0 1];

tspan = [0:0.05:3.95];

[t, y] = ode45(@(t, y) odefcn(t, y), tspan, y0);

function dydt = odefcn(t, y)
    dydt = zeros(7, 1);
    dydt(2) = -0.75 * y(1) * y(3);
    dydt(1) = 0.75 * y(2) * y(3) + 0.2;
    dydt(3) = 0;
    dydt(4) = 1 / 2 * (y(1) * y(7) + y(2) * y(6) - y(3) * y(5));
    dydt(5) = 1 / 2 * (y(2) * y(7) + y(3) * y(4) - y(1) * y(6));
    dydt(6) = 1 / 2 * (y(3) * y(7) + y(1) * y(5) - y(2) * y(5));
    dydt(7) = -1 / 2 * (y(1) * y(4) + y(2) * y(5) + y(3) * y(6));
end

Note that you need dydt = zeros(7, 1) because otherwise the code creates a column vector. Without it, you may run into problems.

ofey
  • 135
  • 8
  • inital_val is just the I.C for this system, but when I run the code it looks right. Thanks a lot!!! I am going to study your code like a religion!! – SHiSA Mar 11 '22 at 20:18
  • One could make this a little closer to the model by using `w=y(1:3); q=y(4:7);` to separate the Euler angles from the orientation quaternion. – Lutz Lehmann Mar 12 '22 at 07:19
2

There are two fundamental problems with your code. First problem is you don't have the correct signature of the derivative function needed for use with ode45( ). Second problem is that blindly integrating quaternion elements this way will lead to non-unit quaternions. I don't know of an easy way to enforce this restriction when using ode45( ).

Let's tackle the easy problem first, the derivative function signature. Also your use of variables inside this function isn't correct. It needs to be like this:

function dy = soe(t,y) % Required signature is (time,state_vector)
% Pick off the named variables from the state vector y
omega1 = y(1);
omega2 = y(2);
omega3 = y(3);
q1 = y(4);
q2 = y(5);
q3 = y(6);
q4 = y(7);
% Calculate the variable derivatives
omega1_dot = -0.75*omega2*omega3; % IS THIS A TYPO? Should it be omega1*omega3?
omega2_dot = 0.75*omega2*omega3+0.2;
omega3_dot = 0;
q1_dot = (1/2)*(q4*omega1+omega2*q3-omega3*q2);
q2_dot = (1/2)*(q4*omega2+omega3*q1-omega1*q3);
q3_dot = (1/2)*(q4*omega3+omega1*q2-omega2*q2);
q4_dot = (1/2)*(q1*omega1+q2*omega2+q3*omega3);
% Group variable derivatives into a single column vector
dy = [omega1dot; omega2dot; omega3dot; q1_dot; q2_dot; q3_dot; q4_dot];
end

Also your angular rate initial values seem completely unreasonable for a rotation problem. They will be interpreted as [1 -1 2] radians/second. I suspect you wanted degrees/second for this, so you will need to scale down those values.

The second problem, that of keeping the quaternion elements forming a unit quaternion, is not going to be easy to overcome using ode45( ). All of the little integration errors will creep in driving the quaternion magnitude away from 1. Any subsequent use of this quaternion in your other code will have problems because of this. My only advice in this is to maybe abandon ode45( ) altogether and write your own custom integrator (e.g., RK4) so that at each step you can re-normalize the quaternion elements.

James Tursa
  • 2,242
  • 8
  • 9
  • Thanks a lot for the help. I definitely see now what I was doing wrong (facepalming atm) and how it works. As for the other things, it was just an ode45 problem, so the values don't mean much. Again thank you so much!! – SHiSA Mar 11 '22 at 20:49
  • The initial values seem not that problematic. 2rad/sec gives a period of about 3, with a time step of 0.05 this are 60 samples per period. This is quite sufficient. – Lutz Lehmann Mar 12 '22 at 07:16