时间戳

首页 > 解决方案 > matlab中点规则

matlab - matlab中点规则

问题描述

您好，我被要求为中点规则创建一个 matlab 代码。我所拥有的是欧拉方法的代码，所以我必须进行一些修改，但我正在努力做到这一点我有以下

function H = heun(f,a,b,ya,M)
h = (b-a)/M;
T = zeros(1,M+1);
Y = zeros(1,M+1);
T = a:h:b;
Y(1) = ya;
for j = 1 : M
k1=feval(f,T(j),Y(j));
k2=feval(f,T(j+1),Y(j)+h*k1);
Y(j+1)=Y(j)+(h/2)*(k1+k2);
end
H = [T' Y'];
function f = dif1(t,y)
f=(t-y)/2;

所以是的，我的函数是 y'=(ty)/2 以及我在命令窗口中定义的间隔和其他内容..

如果可以将 for 循环变成中点规则，我认为这是要走的路，任何帮助将不胜感激。

标签： matlabnumerical-methodsdifferential-equationsrunge-kutta

解决方案

下面是我在 MATLAB 中完成的欧拉方法的一个实现，用于求解一对耦合的一阶 DE。它解决了由以下表示的谐振子：

y1(t+h) = y1(t) + h*y2(t)

y2(t+h) = y2(t) + h*(-A/M y1(t) -B/M y1(t)/|y1(t)|)

% Do the integration using the Euler Method
    while(T<=T1)
        % Update the position of the pt mass using current velocity
        Y1(i+1) = Y1(i) + H*Y2(i);

        % Update the velocity of the pt mass checking if we are turning 
        % within the C/A band
        if ( (abs(Y2(i) < th) && abs(Y1(i)) < C/A) )
            Y2(i+1) = Y2(i);
        else
            Y2(i+1) = Y2(i) + H * ( ((-A/M)*Y1(i)) - (B/M)*sign(Y2(i)) );
        end

        % Incriment the time by H 
        T = T + H;
        % Increase the looping index variable
        i = i + 1;
    end

在没有明确解决您的作业问题的情况下，我希望这个示例有所帮助。需要注意的几件事：在这个特定的例子中，if 语句考虑了静摩擦——所以忽略这一点，只看“else”之后发生的事情。

另请注意，我将初始条件 y1(0) 和 y2(0) 定义为 Y1(i=1) 和 Y2(i=1)，因此从 Yj(i+1) 开始给出 Yj(2)。请注意，T1 是模拟的结束时间。

下面是使用改进欧拉方法的相同问题。如果你推导出这个系统的更新方程，你应该能够浏览代码。

% Do the integration using the Improved Euler Method
    % Ki_j = K^i_j
    while(T<=T1)
        % Calculate K^i_j's

        K1_1 = Y2(i);

        % Must check if we are turning within C/A
        if ( (abs(Y2(i) < th) && abs(Y1(i)) < C/A) )
            K1_2 = 0;
        else
            K1_2 = (-A/M)*Y1(i) - (B/M)*sign(Y2(i));
        end

        K2_1 = Y2(i)+H*K1_2;

        % Checking if we are turning within C/A
        if ( (abs(Y2(i) < th) && abs(Y1(i)) < C/A) )
            K2_2 = 0;
        else
            K2_2 = (-A/M)*(Y1(i) + H*K1_1) - (B/M)*sign(Y2(i)+ H*K1_2);
        end


        % Update the position and velocity
        Y1(i+1) = Y1(i)+ (H/2)*(K1_1+K2_1);
        Y2(i+1) = Y2(i) + (H/2)*(K1_2+K2_2);

        % Incriment the time by H 
        T = T + H;
        % Increase the looping index variable
        i = i + 1;
    end

推荐阅读