правило средней точки для Matlab - PullRequest
Новая
       68

правило средней точки для Matlab

0 голосов
/ 04 октября 2018

Здравствуйте, меня попросили создать код 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 правилом средней точки, я думаю, что это путь, любая помощь будет оценена.

1 Ответ

0 голосов
/ 04 октября 2018

Ниже приведена реализация в 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
...