计算科学 - 使用事件函数在两组 ODE 之间切换 - 吾爱随笔录

我有一堆 ODE，我试图在 MATLAB 中使用 ode45 来解决。我隐藏了方程式的细节以使其简单（以便构建通用算法）！

$\frac{d R}{d t} = F_{1} (R, N) ........................(Eq.1)$ $\frac{dR}{dt}= F_{1}(R,N) \quad\text{........................(Eq.1)}$ $\frac{d N_{1}}{d t} = F_{2} (t) .............................(Eq.2)$ $\frac{dN_{1}}{dt}= F_{2}(t) \quad\text{.............................(Eq.2)}$ $\frac{d N_{2}}{d t} = F_{3} (R, N) .......................(Eq.3)$ $\frac{dN_{2}}{dt}= F_{3}(R,N) \quad\text{.......................(Eq.3)}$ $\frac{d N}{d t} = \frac{d N_{1}}{d t} or \frac{d N_{2}}{d t} based on a criteria of m a x (- \frac{d N_{1}}{d t}, \frac{d N_{2}}{d t})$ $\frac{dN}{dt}= \frac{dN_{1}}{dt} \quad \text{or} \quad \frac{dN_{2}}{dt} \quad \text{based on a criteria of } max(-\frac{dN_{1}}{dt},\frac{dN_{2}}{dt})$ 所以，基本上我想要的是，解决 2 种情况下的 R，N。ODE 求解器应在满足条件时停止，并应通过求解不同的 ODE 集以案例 2 重新启动。 ~~$Case 1 : \frac{d R}{d t}; \frac{d N_{1}}{d t}$ $\text{Case 1}: \frac{dR}{dt}; \frac{dN_{1}}{dt}$ $Case 2 : \frac{d R}{d t}; \frac{d N_{2}}{d t}$ $\text{Case 2}: \frac{dR}{dt}; \frac{dN_{2}}{dt}$~~

我知道我需要一个事件函数来做到这一点。不幸的是，我是新手，请您帮忙（检查我的代码是否正常运行）。我也想不出任何测试函数来测试算法。让我知道是否需要任何其他信息。

clc
clear

y0   = [0.1,0.1,0.1,0.21]; %random initial nos.
t0   = 0;                  %random no.
tEnd = 10;                 %random no.

opt = odeset('RelTol', 1e-3);

y = y0;
t = t0;

State = 1;

while t(end) < tEnd
    yprime = myFunction(t, y, State);
    
    if yprime(2) > yprime(3)
        State = 1;
    else
        State = 2;
    end
    
    
    fcn = @(t,y) myFunction(t, y, State);
    
    opt = odeset('Events', @(t, y) myEvents(t,y,State));
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :),opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end

plot(t,y(:,4))


function [yprime] = myFunction(t, y, State)
yprime = [5*y(1); 4*y(4) + 77*t ; 30*y(4) - 7*t]; %random funcs.
if State == 1
    yprime(4) = yprime(2);
    disp('state is 1')
else
    yprime(4) = yprime(3);
    disp('state is 2')
end
end


function [val,isterm,dir] = myEvents(t,y,State)

fun = myFunction(t, y, State);

val      =  fun(2) - fun(3);
isterm   =  0;
dir      =  0;

end