ODE solver using Lobatto IIIA with table of coefficients

Question

So I'm trying to figure out how to solve a given equation y'=-y+2te^(-t+2) for t in [0,10], step of 0.01 and y(0)=0.

I am supposed to solve it using the Lobatto IIIA method following a Butcher tableau:

Coefficients table

So far, this is what I got:

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

It obviously makes no sense from the point when I equal f2 to itself, when the variable is still undefined.

Any help would be much appreciated :)

Cheers

You will need a predictor-corrector loop or apply Newton's or a similar method to solve the implicit equations. — Lutz Lehmann, Jun 05 '18 at 04:58

score 0 · Answer 1 · answered Jun 14 '18 at 15:43

You will need a predictor-corrector loop, in the simple case the corrector uses the slope-equations as basis of a fixed-point iteration. In the code below I use also the value of an explicit Euler step, in principle you could initialize all slopes with 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

The graph shows that the result (blue) is qualitatively correct, the exact solution curve (green) is shifted so that two distinct curves can be seen.

ODE solver using Lobatto IIIA with table of coefficients

1 Answers1