1

在 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 新手，非常欢迎任何帮助！

4

2 回答 2

3

首先，您应该使用subs来评估一个符号函数。另一种方法是使用matlabFunction将所有符号表达式转换为匿名函数，正如 Horchler 所建议的那样。

其次，您正在集成ODE，就好像它是dx/dt. 如果您对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 双精度浮点格式可以毫不费力地表示2.225073858507201e-308和1.797693134862316e+308(realmin和realmax) 之间的数字。因此，对于您给出的系数 (O(10 14 ))，这绝对不是问题。如果您不采取预防措施（重新调整为[-1 +1]，以不同的单位重新表述问题，...），您可能会损失几位数的精度，但与的算法错误ode45。

<RANDOM_OPINIONATED_RANT>

第四，你为什么要使用符号数学来达到这个目的？！您正在进行数值积分，这意味着无论如何都没有解析解。那为什么还要为符号而烦恼呢？与符号进行集成（通过vpa偶数）将是几十，几百，是的，通常甚至比保持（或重新实现）所有数字（有些人认为在 MATLAB 中与裸机相比已经很慢）慢数千倍。金属方法）。

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

  • 使用符号进行推导、证明定理、简化表达式……
  • 使用数字来实现任何期望实际数字的算法或函数。

换句话说，用于绘图的符号，用于处理的数字。并且恰好零符号应该出现在任何算法的任何良好实现中。

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

对我来说，它与邪恶有相似之处，为什么它应该eval有类似的原因。是。避免。

</RANDOM_OPINIONATED_RANT>

于 2017-01-03T15:06:06.687 回答
0

With the full code, it's easy to come up with a complete solution: 

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

Results:

position as a function of time velocity as a function of time

Note that I've not used symbolics anywhere, except to double-check my hand-derived derivatives.

于 2017-01-03T17:09:03.430 回答