用matlab求解单循环的微分方程-Java 学习之路

我的机械系统有以下等式：

xdot = Ax+ Bu

我想在一个循环中解决这个等式，因为在每个步骤中我都需要更新你，但求解器如 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

我是否有另一种方法可以在一个循环中解决这个等式一次（非单一时间步） .

1 回答

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

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

您会注意到我使用的是逻辑索引而不是 if -statement . if -statements中的向量仅在所有元素都为真时才为真，反之亦然为假 . 因此，我创建一个逻辑数组并使用它来分母 1 或 Alpha ，具体取决于 x / U 的每一行的信号值 .

'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

回复于 2024-04-20T17:48:45+08:00

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

1 回答

相关问题