ODE45 以非常大的数字作为约束

Question

2^nd 要在 MATLAB 中求解的 ODE：

  ( (a + f(t))·d²x/dt² + (b/2 + k(t))·dx/dt ) · dx/dt - g(t) = 0

边界条件：

 dx/dt(0) = v0

哪里

t是时间，
x是位置
dx/dt是速度
d2x/dt2是加速度
a、b、v0是常量
f(t)、k(t) 和 h(t) 是依赖于 t 的已知函数（我不写了，因为它们很大）

举个例子，使用符号变量：

syms t y
%% --- Initial conditions ---
   phi = 12.5e-3;
   v0 = 300;            
   e  = 3e-3;           
   ro = 1580;          
   E  = 43e9;   
   e_r = 0.01466;  
   B = 0.28e-3; 

%% --- Intermediate calculations ---
   v_T = sqrt(((1 + e_r) * 620e6) /E) - sqrt(E/ro) * e_r;
   R_T = v_T * t;                   
   m_acc = pi * e * ro *(R_T^2);  
   v_L = sqrt (E/ro);           
   R_L = v_L * t;    
   z = 2 * R_L;                   
   E_4 = B * ((e_r^2)* B * (0.9^(z/B)-1)) /(log(0.9));
   E_1 = E * e * pi * e_r^2 * (-phi* (phi - 2*v_T*t)) /16; 
   E_2 = pi * R_T^2 * 10e9;        
   E_3 = pi * R_T^2 * 1e6 * e;   

    %% Resolution of the problem  
        g_t = -diff(E_1 + E_2 + E_3, t); 

f(t,y)=(g_t - (pi*v_T*e*ro/2 + E_4) * y^2 /(y* (8.33e-3 + m_acc))];
fun=matlabFunction(f);
[T,Y]=ode45(fun,[0 1], v0]);

我如何重写它以得到 x 为 y=dx/dt？我是 Matlab 的新手，非常欢迎任何帮助！

Answer 1

首先，你应该使用subs to evaluate a symbolic function. Another approach is to use matlabFunction将所有符号表达式转换为匿名函数，正如Horchler所建议的那样。

Second，您正在对 ODE 进行积分，就好像它是 dx/dt 中的 1^st 阶一样。如果您对 x(t) 和 dx/dt(t) 感兴趣，那么您必须像这样修改函数：

fun = @(t,y) [y(2);   
             ( subs(g) - (b/2 + subs(k))*y(2)*y(2) ) / ( y(2) * (a + subs(f))) ];

当然，还要为 x0 = x(0) 和 v0 = dx/dt(0) 提供初始值。

第三，参数的绝对值几乎从来都不是真正关心的问题。 IEEE754 double-precision floating point format 可以毫不费力地表示 2.225073858507201e-308 和 1.797693134862316e+308 之间的数字（分别为 realmin 和 realmax）。所以对于你给的系数（O(10¹⁴)），这绝对不是问题。如果您不采取预防措施（重新缩放到 [-1 +1]，以不同的单位重新表述问题，...），您可能会失去几位数的精度，但由此产生的相对误差很可能很小与 ode45 所犯的算法错误相比微不足道。

<RANDOM_OPINIONATED_RANT>

第四,为什么你为了这个目的使用符号数学？！您正在进行 numerical 积分，也就是说，无论如何都没有解析解。那为什么还要用符号呢？与符号进行集成（通过 vpa 甚至）将比保持（或重新实现）所有数值慢数十、数百，是的，通常甚至千倍（有些人认为与裸机方法相比，MATLAB 中的速度已经很慢了）。

是的，当然，对于这个特定的、单独的、孤立的用例，它可能并不重要，但为了将来，我强烈建议您学习：

使用符号进行推导、证明定理、简化表达式...
使用数值实现任何需要实际数字的算法或函数。

换句话说，symbolics for drafting，numerics for 嘎吱作响。并且 zero 符号应该出现在任何算法的任何良好实现中。

虽然可能在某种程度上混合使用它们，但这并不意味着这样做是个好主意。事实上，这几乎从来没有。它是唯一可行选择的少数孤立案例并不能证明这种方法是正确的。毕竟，它们是罕见的孤立案例，远非丰富的常态。

对我来说它与邪恶有相似之处eval, with similar reasons for why it Should. Avoided。

</RANDOM_OPINIONATED_RANT>

Answer 2

有了完整的代码，很容易想出一个完整的解决方案：

% Initial conditions
phi = 12.5e-3;
v0  = 300;
x0  = 0;    % (my assumption)
e   = 3e-3;
ro  = 1580;
E   = 43e9;
e_r = 0.01466;
B   = 0.28e-3;

% Intermediate calculations
v_T = sqrt(((1 + e_r) * 620e6) /E) - sqrt(E/ro) * e_r;
R_T = @(t) v_T * t;
m_acc = @(t) pi * e * ro *(R_T(t)^2);
v_L = sqrt (E/ro);
R_L = @(t) v_L * t;
z   = @(t) 2 * R_L(t);
E_4 = @(t) B * ((e_r^2)* B * (0.9^(z(t)/B)-1)) /(log(0.9));

% UNUSED
%{
E_1 = @(t) -phi * E * e * pi * e_r^2 * (phi - 2*v_T*t) /16;
E_2 = @(t) pi * R_T(t)^2 * 10e9;
E_3 = @(t) pi * R_T(t)^2 * 1e6 * e;
%}

% Resolution of the problem
g_t = @(t)  -( phi * E * e * pi * e_r^2 * v_T / 8 + ...  % dE_1/dt
               pi * 10e9 * 2 * R_T(t) * v_T + ...        % dE_2/dt
               pi * 1e6 * e * 2 * R_T(t) * v_T );        % dE_3/dt

% The derivative of Z = [x(t); x'(t)] equals Z' = [x'(t);  x''(t)]
f = @(t,y)[y(2);
          (g_t(t) - (0.5*pi*v_T*e*ro + E_4(t)) * y(2)^2) /(y(2) * (8.33e-3 + m_acc(t)))];

% Which is readily integrated
[T,Y] = ode45(f, [0 1], [x0 v0]);

% Plot solutions
figure(1)
plot(T, Y(:,1)) 
xlabel('t [s]'), ylabel('position [m]')

figure(2)
plot(T, Y(:,2)) 
xlabel('t [s]'), ylabel('velocity [m/s]')

结果：

请注意，我没有在任何地方使用符号，除非仔细检查我的手推导数。

ODE45 以非常大的数字作为约束

ODE45 with very large numbers as constraints

matlab

symbolic-math

ode

ode45