2

我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<3600if

在此处输入图像描述

4

3 回答 3

2

这是一种常见的问题。用户经常希望不连续地更改 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'根据上一次调用中使用的最后一步大小指定一个选项来帮助求解器。通过在命令窗口中键入来查看ballode演示以获取更多详细信息。edit ballode

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

于 2013-12-06T17:09:31.770 回答
1

如果你有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到介于0 和 1 之间Q1，Q2而不是 0 和 1。

于 2013-12-06T17:08:42.463 回答
1

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

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

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

于 2013-12-06T17:32:36.377 回答