不带计数器的自定义周期函数答案

【问题标题】：Custom periodic function without counter不带计数器的自定义周期函数
【发布时间】：2013-12-24 02:37:19
【问题描述】：

我正在使用ode45 解决一个简单的ODE：

function dCdt=u_vent(t,C)
if t>   600 &&  t<= 720
    Q=Q2;               
elseif t>   1320    &&  t<= 1440
    Q=Q2;           
elseif t>   2040    &&  t<= 2160
    Q=Q2;           
elseif t>   2760    &&  t<= 2880
    Q=Q2;               
elseif t>   3480    &&  t<= 3600
    Q=Q2;
else
    Q=Q1;
end

V=100;
C_i=400;
S=100;

dCdt=Q/V*C_i+S/V-Q/V*C(1);
return

我用then来解决：

[t,C]=ode45(@u_vent, [0 1*3600], 400);

我想为Q、0<t<3600 创建一个周期函数，例如图片中的那个，而不使用那些if 语句...有什么想法吗？

【问题讨论】：

如果你有 Simulink，那么它很容易以图形方式完成
我假设您知道您的系统有解析解？如果你不擅长微分方程，你可以使用Matlab的dsolve函数：syms Q V C_i S t C(t) C0;dsolve(diff(C,t)==(Q*(C_i-C)+S)/V,C(0)==C0)。

标签： matlab differential-equations numerical-integration

【解决方案1】：

这是一种常见的问题。用户经常希望不连续地更改 Matlab 的 ODE 求解器使用的积分函数中的变量。不幸的是，这通常是个坏主意。这些求解器假设微分方程是平滑且连续的。充其量你的代码会更慢。在最坏的情况下，你会有更大的错误甚至完全错误的结果。使用诸如ode15s 之类的刚性求解器可能会有所帮助，但它也假定了连续性。由于您指定的系统具有解析解，因此“模拟”它的最简单方法实际上可能是通过这个时间函数传递脉冲序列。

参数在时间上不连续变化的情况，即相对于自变量，很容易解决。只需在每个时间跨度上整合您想要的时段数：

t1 = 600;
t2 = 120;
TSpan = [0 t1]; % Initial time vector
NPeriods = 5;   % Number of periods
C0 = 400;       % Initial condition
Q1 = ...
Q2 = ...
V = 100;
C_i = 400;
S = 100;
u_vent = @(t,C,Q)(Q*(C_i-C(1))+S)/V; % Integration function as anonymous function
Cout = C0;       % Array to store solution states
tout = Tspan(1); % Array to store solution times
for i = 1:NPeriods
    [t,C] = ode45(@(t,C)u_vent(t,C,Q1), TSpan, C0);
    tout = [tout;t(2:end)];   % Append time, 2:end used to avoid avoid duplicates
    Cout = [Cout;C(2:end,:)]; % Append state
    TSpan = [0 t2]+t(end);    % New time vector
    C0 = C(end);              % New initial conditions
    [t,C] = ode45(@(t,C)u_vent(t,C,Q2), TSpan, C0);
    tout = [tout;t(2:end)];
    Cout = [Cout;C(2:end,:)];
    TSpan = [0 t1]+t(end);
    C0 = C(end);
end

这允许ode45 根据一组初始条件在指定时间内积分 ODE。然后在上一次积分结束时重新开始积分，新的不连续初始条件或不同的参数。这导致瞬态。您可能不喜欢代码的外观，但这就是它的完成方式。如果参数相对于状态变量不连续地变化，那就更棘手了。

可选。每次您调用/重新启动ode45 时，该函数都必须确定要使用的初始步长。这是重新启动求解器的主要（最小）成本/开销。在某些情况下，您可以根据上一次调用中使用的最后一步大小指定 'InitialStep' 选项来帮助求解器。通过在命令窗口中输入edit ballode 来查看ballode 演示以获取更多详细信息。

一个注释。 tout 和 Cout 数组在每个集成步骤之后附加到它们自身。这实际上是动态内存分配。只要NPeriods 相当小，这可能不会成为问题，因为最近版本的 Matlab 中的动态内存分配可能非常快，而且您只需在大块中重新分配几次。如果您指定固定步长输出（例如TSpan = 0:1e-2:600;），那么您将确切知道要为tout 和Cout 预分配多少内存。

【讨论】：

啊，太棒了！非常感谢！这似乎是唯一明智的解决方案。如果 Q 是一个函数，例如 @am304 建议的函数，那么同样的事情是否仍然可行？
如果它是一个平滑变化的输入（例如，正弦强制），您可以将其包含在积分函数中。像您这样的非平滑输入可能有效，也可能无效。 “可能有效”是指得到不完全错误但仍有更多错误的结果。这取决于特定的 ODE 和求解器以及求解器的设置。一些 ODE 在受到瞬变或错误会爆炸时需要仔细集成。无论如何，如果您使用此答案来集成具有不连续性的系统（对于相同的求解器容差），错误将更大。这是一个坏习惯。

【解决方案2】：

如果你有Signal Processing Toolbox，你可以使用类似的东西：

>> t = 0:3600;
>> y = pulstran(t,[660:720:3600],'rectpuls',120);
>> plot(t,y)
>> ylim([-0.1 1.1])

给出以下内容（在 Octave 中，在 MATLAB 中应该相同）：

然后您需要将y 缩放到Q1 和Q2 之间，而不是0 和1。

【讨论】：

horchler 的建议可能比我的好。

【解决方案3】：

不一定是最好的方法，因为连续性假设仍然被打破，但是在没有 if 链的情况下生成矩形脉冲序列的方法是：

Q = Q2 + (Q1 - Q2) * (mod(t, period) < t_rise);

适用于标量和向量上下文。

【讨论】：