ODE решатель с использованием Lobatto IIIA с таблицей коэффициентов - PullRequest
Новая
       71

ODE решатель с использованием Lobatto IIIA с таблицей коэффициентов

0 голосов
/ 04 июня 2018

Так что я пытаюсь выяснить, как решить данное уравнение y '= - y + 2te ^ (- t + 2) для t в [0,10], шаг 0,01 иy (0) = 0.

Я должен решить эту проблему, используя метод Лобатто IIIA по таблице Мясника:

Таблица коэффициентов

Итакдалеко, вот что я получил: 

function lob = Lobatto_solver()
     h = 0.01;
     t = 0:h:10;
     y = zeros(size(t));
     y(1) = 0;

     f = @(t,y) -y + (2*t*(exp(-t+2)));

     % Lobatto IIIA Method
     for i=1:numel(y)-1
         f1 = f(t(i), y(i));
         f2 = f(t(i)+(1/2)*h, y(i) + (5/24)*h*f1 + (1/3)*h*f2 + (-1/24)*h*f3);
         f3 = f(t(i)+h, y(i) + (1/6)*h*f1 + (2/3)*h*f2 + (1/6)*h*f3);
         y(x) = y(i) + h*((-1/2)*f1 + (2)*f2 + (-1/2)*f3);
     end
 end

Это, очевидно, не имеет смысла с того момента, когда я приравниваю f2 к себе, когда переменная все еще не определена.

Любая помощь будет многооценил:)

Ура

1 Ответ

0 голосов
/ 14 июня 2018

Вам потребуется цикл предиктор-корректор, в простом случае корректор использует уравнения наклона в качестве основы итерации с фиксированной точкой.В приведенном ниже коде я также использую значение явного шага Эйлера, в принципе вы можете инициализировать все уклоны с помощью f1. 

function lob = Lobatto_solver()
    h = 0.01;
    t = 0:h:10;
    y = zeros(size(t));
    y(1) = 0;

    f = @(t,y) -y + (2*t*(exp(-t+2)));

    % Lobatto IIIA Method
    for i=1:numel(y)-1
        f1 = f(t(i), y(i));
        f3 = f(t(i)+h, y(i)+h*f1)
        f2 = (f1+f3)/2;
        for k=1:3
          f2 = f(t(i)+(1/2)*h, y(i) + (5/24)*h*f1 + (1/3)*h*f2 + (-1/24)*h*f3);
          f3 = f(t(i)+h, y(i) + (1/6)*h*f1 + (2/3)*h*f2 + (1/6)*h*f3);
        end;
        y(i+1) = y(i) + h*((-1/2)*f1 + (2)*f2 + (-1/2)*f3);
    end
    plot(t,y,t,0.05+t.^2.*exp(-t+2))
end

График показывает, что результат (синий) качественно правильный,кривая точного решения (зеленая) смещена так, что можно увидеть две различные кривые.enter image description here

...