如何提取与 ode45 一起使用的函数中的特定值？

问题描述 投票：0回答：1

我有一个与 ode45 一起使用的函数。函数返回速度和加速度，我用ode45求解得到位置和速度。

在函数内，我通过将加速度乘以质量来定义“力”，我想将其提取到函数外部以绘制它。

我尝试过两种方法。

第一个方法是这个（代码在下面），但是它在矩阵计算中返回错误。

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

矩阵

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)]

像这样在每个 2x3 数组中仅存储第一列会出现问题。

...
   -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

必须得到这样的精确结果。

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

第二种方法涉及使用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);

运行没有错误，但给出了错误的结果。

我用的cellfun有问题吗？任何援助将不胜感激。谢谢你。

我的完整代码

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

matlab cell ode45
1个回答
0
投票

如果您将 

constrained_ECI
与 
ode45
一起使用，那么您需要匹配 
ode45
所需的格式（请参阅文档），并且不要添加额外的输出。解决这个问题的方法很少。

也许最简单、最好的方法就是创建两个单独的函数：一个用于数值积分，另一个用于根据积分结果计算所需的值。这样做的优点是不会在 ODE 函数中包含额外的不需要的计算。如果您确实需要 ODE 中另一个函数的结果，您可以根据需要使用第二个函数。这是一个使用更简单的 ODE 系统的示例：

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

另一种选择是在集成函数中包含一些逻辑来切换它返回的内容。但请注意，这可能会在求解 ODE 时增加一些额外的计算开销。通常，人们不希望在 ODE 函数中包含任何分支，但这里在积分过程中只应调用其中一个分支。例如：

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

请注意，在上述两种情况下，如果您要传递来自 

ode45
的完整结果向量，则需要对函数进行一些向量化，因为名义上 ODE 积分函数仅期望缩放时间和状态向量那个时候。

最后，有一些可以使用

global
变量的方案可能有效，但由于性能下降和糟糕的编码风格，通常不推荐。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.