Matlab中ODE直到平衡的求解系统

Question

function [t,N] = SI1(s0,i0,Tmax,b,d,<more parameters>)

% Initial conditions
y0 = [s0 i0];

% parameters
P(1) = b;      
P(2) = d;      
P(3) = ... etc

% time steps
time = [0 Tmax];

% solve ODE system
[t,N] = ode45(@(t,Y) SI1ode(t,Y,P), time, y0);

end


function Y = SI1ode(t,x,P)

%parameters
b = P(1);      
d = P(2);      
...etc

%State variables
S = x(1);
I = x(2);

Y(1) = equation 1;
Y(2) = equation 2;
Y=Y';   %make column vector
end

而我的行书是：

Tmax = 100;
s0 = 100;
i0 = 1e-2;

b = 2;      
...etc

D = [0 .01 0.05 0.1 0.15 0.2 0.3 0.4 0.5 .6 .7 .8 .9 1];

for j=1:length(D)

    d = D(j);
    [t,N] = SI1(s0,i0,Tmax,b,d,<more parameters>);
    result(j) = N(end,2)/(N(end,1)+N(end,2));

end

Answer 1

可以通过使用输出函数或事件来中断集成。如果您更喜欢一个事件，您应该编写类似于If equality,isterminal(1)=1；end，而对于输出函数，如果平衡，status=1，则代码将变为；结束。

关于如何测试平衡条件，您应该将解的前一个值与当前值进行比较，这是一个可以逐个组件进行的操作，也可以使用规范标准进行操作。在这两种情况下，这都类似于用有限差分计算解的导数。

此操作意味着您必须有权访问在前面的时间步骤中计算出的解决方案。可以使用output函数将这些值保存在全局变量中。