用matlab循环求解单次微分方程

Question

我有一个具有以下方程式的机械系统：

xdot = Ax+ Bu

我想在一个循环中求解这个方程，因为在每一步中我都需要更新 u，但是像 ode45 或 lsim 这样的求解器会在一个时间间隔内求解微分方程。

 for i = 1:10001
    if  x(i,:)>= Sin1 &  x(i,:)<=Sout2
        U(i,:) = Ueq - (K*(S/Alpha))
    else
        U(i,:) = Ueq - (K*S)
    end
   % [y(i,:),t,x(i+1,:)]=lsim(sys,U(i,:),(time=i/1000),x(i,:));
   or %[t,x] = ode45(@(t,x)furuta(t,x,A,B,U),(time=i/1000),x)
end

我是否有另一种方法可以在一个循环中一次性求解这个方程（不是单个时间步）。

Answer 1

有许多方法可以跨函数调用更新和存储数据。 对于 ODE 套件，我开始喜欢所谓的“闭包”。 闭包基本上是一个嵌套函数，从其父函数访问或修改变量。

下面的代码通过将传递给ode45 的右侧函数和'OutputFcn' 包装在名为odeClosure() 的父函数中来利用此功能。

您会注意到我使用的是逻辑索引而不是if 语句。 if-statements 中的向量只有在所有元素都为真时才为真，反之亦然。 因此，我创建了一个逻辑数组，并根据x/U 的每一行的信号值，使用它来使分母为1 或Alpha。

'OutputFcn'storeU() 在成功的时间步之后由ode45 调用。 该函数增长U 存储阵列并适当地更新它。 数组U 将具有与tspan 请求的解点数相同的列数（在这个虚构的示例中为12）。 如果一个成功的整步跳过任何请求的点，则调用该函数，其中包含所有请求的中间时间及其相关的解值（因此x 可能是矩形，而不仅仅是一个向量）；这就是为什么我在storeU 中使用bsxfun 而不是rhs。

示例函数：

function [sol,U] = odeClosure()

    % Initilize
%     N     = 10          ;
    A     = [ 0,0,1.0000,0; 0,0,0,1.0000;0,1.3975,-3.7330,-0.0010;0,21.0605,-6.4748,-0.0149];
    B     = [0;0;0.6199;1.0752 ] ;
    x0    = [11;11;0;0];
    K     = 100;
    S     = [-0.2930;4.5262;-0.5085;1.2232];
    Alpha = 0.2          ;
    Ueq   = [0;-25.0509;6.3149;-4.5085];
    U     = Ueq;
    Sin1  = [-0.0172;-4.0974;-0.0517;-0.2993];
    Sout2 = [0.0172 ; 4.0974; 0.0517; 0.2993];

    % Solve
    options = odeset('OutputFcn', @(t,x,flag) storeU(t,x,flag));
    sol     = ode45(@(t,x) rhs(t,x),[0,0.01:0.01:0.10,5],x0,options);


    function xdot = rhs(~,x)

        between = (x >= Sin1) &  (x <= Sout2);
        uwork   = Ueq - K*S./(1 + (Alpha-1).*between);
        xdot    = A*x + B.*uwork;

    end

    function status = storeU(t,x,flag)

        if isempty(flag)
            % grow array
            nAdd      = length(t)           ;
            iCol      = size(U,2) + (1:nAdd);
            U(:,iCol) = 0                   ;

            % update U
            between   = bsxfun(@ge,x,Sin1) & bsxfun(@le,x,Sout2);
            U(:,iCol) = Ueq(:,ones(1,nAdd)) - K*S./(1 + (Alpha-1).*between);
        end

        status = 0;
    end

end