matlabcellode45

How to extract specific values in a function used with ode45?

I have a function for use with ode45. The function returns velocity and acceleration, and I use ode45 to solve it to obtain position and velocity.

within the function, I define "Force" by multiplying acceleration with mass, and I want to extract this, outside the function to plot it.

I have tried two methods.

The first method is this (code down below), but it returns errors in matrix calculations.

[~, Force] = constrained_ECI(t, result.', G, M, i, ohm, n)

matrix

R_bar = [- sin(t.*nval).*cos(ohmval) - cos(t.*nval).*cos(ival).*sin(ohmval), cos(t.*nval).*cos(ival).*cos(ohmval) - sin(t.*nval).*sin(ohmval), cos(t.*nval).*sin(ival); ...
        sin(ival)*sin(ohmval), -cos(ohmval)*sin(ival), cos(ival)]

an issue with storing only the first columns in each 2x3 array like this.

...
   -0.0846    0.1517    0.9848
   -0.0854    0.1512    0.9848
   -0.0858    0.1510    0.9848
   -0.0861    0.1508    0.9848
   -0.0865    0.1506    0.9848
   -0.0868    0.1504    0.9848
    0.4924   -0.8529    0.1736

have to get exact result like this.

...
   -0.2607   -0.3349   -0.9055
    0.4924   -0.8529    0.1736
   -0.2607   -0.3349   -0.9055
    0.4924   -0.8529    0.1736
...

The second method involves using cellfun,

[~,Force] = cellfun(@(t,x) constrained_ECI(t, result.', G, M, i, ohm, n), num2cell(t), mat2cell(result, ones(size(result, 1), 1), size(result, 2)), 'UniformOutput', false)
ForceData = cell2mat(Force);

which runs without errors but gives wrong results.

Is there a problem with the cellfun I used? Any assistance would be appreciated. Thank you.

My Full Code

clc; clear;

%% parameters
G = 6.67430e-11;  % J*m/kg^2 Gravitational constant
M = 5.9722e24; % kg, Earth Mass
rho = 50000;
w0 = 0;
i = deg2rad(80);
ohm = deg2rad(30);
r0 = 7e6; % radius of chief
r0vector = [r0;0;0];

%% time
n = 0.001078; % rad/s
p = 2*pi/n; % period
tspan = [0 p*3];  % Time span for the simulation
opts = odeset('RelTol',1e-12,'AbsTol',1e-14); 

%% Initial conditions
% x0 y0 z0 xdot0 ydot0 zdot0
alpha0 = 0.55;
x0 = rho*sin(alpha0)/2;                                       
y0 = rho*cos(alpha0);                                               
z0 = rho*sin(alpha0);  
xdot0 = rho*n*cos(alpha0)/2;                                          
ydot0 = -rho*n*sin(alpha0);  
zdot0 = rho*n*cos(alpha0); 

xyz0 = [x0;y0;z0];
xyzdot0 = [xdot0;ydot0;zdot0];

hillxyz0 = xyz0 + r0vector;

%% rotation matrix
syms nval time ival ohmval real
t0 = 0;
wval = w0 + nval*(time-t0);

R3w = [cos(wval) sin(wval) 0; -sin(wval) cos(wval) 0; 0 0 1];
R1i = [1 0 0; 0 cos(ival) sin(ival); 0 -sin(ival) cos(ival)];
R3ohm = [cos(ohmval) sin(ohmval) 0; -sin(ohmval) cos(ohmval) 0; 0 0 1];

R_cal = R3w * (R1i * R3ohm);

% hill x y z -> ECI X Y Z
R_0 = subs(R_cal, [nval, ival, ohmval, time], [n, i, ohm, 0]);
ECIXYZ0 = inv(R_0)*hillxyz0;

% hill x' y' z' -> ECI X' Y' Z'
Rdot = diff(R_cal, time);
Rdot0 = subs(Rdot, [nval, ival, ohmval, time], [n, i, ohm, 0]);
ECIXYZdot0 = inv(R_0) * (xyzdot0 - Rdot0 * ECIXYZ0);

initial_state = double(vpa([ECIXYZ0; ECIXYZdot0], 4));

%% ===== Solve =====
[t, result] = ode45(@(t, y) constrained_ECI(t, y, G, M, i, ohm, n), tspan, initial_state, opts);

X = result(:, 1);
Y = result(:, 2);
Z = result(:, 3);
ECImatrix = [X Y Z];

%% ===== method 1 =====
[~, Force] = constrained_ECI(t, result.', G, M, i, ohm, n)

%% ===== method 2 =====
[~,Force] = cellfun(@(t,x) constrained_ECI(t, result.', G, M, i, ohm, n), num2cell(t), mat2cell(result, ones(size(result, 1), 1), size(result, 2)), 'UniformOutput', false)
ForceData = cell2mat(Force);

%% ===== function =====
function [state, Force] = constrained_ECI(t, y, G, M, ival, ohmval, nval)

    m = 1;

    vel = [y(4);y(5);y(6)]; % X dot
    r = [y(1);y(2);y(3)]; % X

    R_bar       = zeros(2,3);
    R_bar_dot   = zeros(2,3);
    R_bar_ddot  = zeros(2,3);
    R_tilde     = zeros(2,3);
    R_tilde_dot = zeros(2,3); 
    R_tilde_ddot= zeros(2,3);

    R_bar = [- sin(t*nval)*cos(ohmval) - cos(t*nval)*cos(ival)*sin(ohmval), cos(t*nval)*cos(ival)*cos(ohmval) - sin(t*nval)*sin(ohmval), cos(t*nval)*sin(ival); ...
        sin(ival)*sin(ohmval), -cos(ohmval)*sin(ival), cos(ival)]
    R_bar_dot = [-nval*cos(ohmval)*cos(t*nval)+nval*sin(ohmval)*cos(ival)*sin(t*nval), -nval*cos(ohmval)*cos(ival)*sin((t*nval))-nval*sin(ohmval)*cos((t*nval)), -nval*sin(ival)*sin(t*nval); ...
        0 0 0];
    R_bar_ddot = [(nval^2)*cos(ohmval)*sin(t*nval)+(nval^2)*sin(ohmval)*cos(ival)*cos(t*nval), -(nval^2)*cos(ohmval)*cos(ival)*cos(t*nval)+(nval^2)*sin(ohmval)*sin(t*nval), -(nval^2)*sin(ival)*cos(t*nval); ...
        0 0 0];
    
    R_tilde = [cos(t*nval)*cos(ohmval) - sin(t*nval)*cos(ival)*sin(ohmval), cos(t*nval)*sin(ohmval) + sin(t*nval)*cos(ival)*cos(ohmval), sin(t*nval)*sin(ival); ...
        sin(ival)*sin(ohmval), -cos(ohmval)*sin(ival), cos(ival)];
    R_tilde_dot = [-nval*cos(ohmval)*sin(t*nval)-nval*sin(ohmval)*cos(ival)*cos(t*nval), -nval*sin(ohmval)*sin(t*nval)+nval*cos(ohmval)*cos(ival)*cos(t*nval), nval*sin(ival)*cos(t*nval); ...
        0 0 0];
    R_tilde_ddot = [-(nval^2)*cos(ohmval)*cos(t*nval)+(nval^2)*sin(ohmval)*cos(ival)*sin(t*nval), -(nval^2)*sin(ohmval)*cos(t*nval)-(nval^2)*cos(ohmval)*cos(ival)*sin(t*nval), -(nval^2)*sin(ival)*sin(t*nval); ...
        0 0 0];

    %% [A11 A12 A13]
    Arow1 = (r')*((R_bar')*R_bar);
    A11 = Arow1(1);
    A12 = Arow1(2);
    A13 = Arow1(3);

    %% [A21 A22 A23]
    Arow2 = [2, -1]*R_tilde;
    A21 = Arow2(1);
    A22 = Arow2(2);
    A23 = Arow2(3);

    %% [b1 b2]
    b1 = -(r')*(R_bar')*R_bar_ddot*(r) - 2*(r')*(R_bar')*R_bar_dot*(vel)...
        -(r')*(R_bar_dot')*R_bar_dot*(r) - 2*(r')*(R_bar_dot')*R_bar*(vel)...
        -(vel')*(R_bar')*R_bar*(vel);
    b2 = -[2 -1]*R_tilde_ddot*(r) - [4 -2]*R_tilde_dot*(vel);

    Amatrix = [A11 A12 A13;A21 A22 A23];
    Bmatrix = [b1;b2];

    A1MPinverse  = (1/(A11^2+A12^2+A13^2))*(Amatrix(1,:)');
    Bbeta = (A11*A21 + A12*A22 + A13*A23)/(A11^2+A12^2+A13^2);
    cMPinverse = (1/(A21^2+A22^2+A23^2-Bbeta^2*(A11^2 + A12^2+A13^2)))*((Amatrix(2,:)')-Bbeta*(Amatrix(1,:)'));

    AMPinverse = [(A1MPinverse - Bbeta*cMPinverse) cMPinverse];

    rnorm = norm(r);

    unconstrained_accel = -(G*M/rnorm^3)*r;
    constrained_accel = (unconstrained_accel) ...
        + (AMPinverse*Bmatrix) ...
        + ((G*M/rnorm^3)*AMPinverse*Amatrix*r);

    %% state, Force
    Force = m*(constrained_accel - unconstrained_accel)
    state = [vel;constrained_accel];

end

correct graph

Solution

  • If you're using constrained_ECI with ode45 then you need to match the format expected by ode45 (see the documentation) and not add extra outputs. There are few ways around this.

    Maybe the simplest and best way to do this is to just create two separate functions: one to use for numerical integration and the other to calculate the value you want base on the integration results. This has the advantage not including additional unneeded computations within your ODE function. If you do need the result of this other function inside your ODE you can the second function there if you like. Here's an example of this with a simpler ODE system:

    c = -0.51;
tspan = [0 12];
y0 = [0; 1; 1];
[t,y] = ode45(@(t,y)f1(t,y,c),tspan,y0);
some_other_result1 = other_result1_func(t,y,c);

function out = f1(t,y,c)
    out = [y(2)*y(3); -y(1)*y(3); other_result1_func(t,y.',c)];
end

function out = other_result1_func(t,y,c)
    out = c*y(:,1).*y(:,2);
end

    Another option is to include some logic within your integration function to switch what it returns. Note however that this may add a bit of extra computation overhead when solving the ODE. Normally, one wouldn't want to include any branching within an ODE function but here only one of the branches should be called during integration. For example:

    c = -0.51;
tspan = [0 12];
y0 = [0; 1; 1];
[t,y] = ode45(@(t,y)f2(t,y,c,false),tspan,y0);
some_other_result2 = f2(t,y.',c,true);

function out = f2(t,y,c,is_other_result)
    other_result = c*y(1,:).*y(2,:)
    if is_intermediate_result
        out = other_result.';
    else
        out = [y(2)*y(3); -y(1)*y(3); other_result];
    end
end

    Note that in both of the above cases, you'll need to do some vectorization of your functions if you with to pass in the full results vectors from ode45 as nominally the ODE integration function only expects a scale time and state vector at that time.

    Finally, there are schemes that could employ global variables that might work but are generally not recommended due to performance hits and poor coding style.