我有一个与 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
如果您将
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
变量的方案可能有效,但由于性能下降和糟糕的编码风格,通常不推荐。