Solving 7 equations of first order equations ode45 Matlab-CodePudding

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

CodePudding user response：

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.

CodePudding user response：

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.