Implicit Trapezoidal Integration Method for ODE

[Namir]

Implicit Trapezoidal Integration Method for ODE
================================================

Given f(x,y) and (x0,y0) as the initial point, xf as the final point, h is the integration step, toler as the tolerance value.

The function to solve for y1 is:

F(y1) = y1 - y0 - h/2*(f(x0,y0) + f(x1,y1))

The loop to calculate y at x=xf is:

Code:

While x0 < xf
  x1 = x0 + h
  y1 = y0 + h*f(x0,y0)
  diff = 2*toler
  i = 1
  max = 100
  While abs(diff) >= toler and i <= max
    i += 1
    hx = 0.01*(1 + abs(y1)
    F0 = F(y1)
    Fp = F(y1+hx)
    Fm - F(y1-hx)
    diff = 2*hx*F0/(Fp - Fm)
    x -= diff
  end
  y0 = y1
  x0 = x1
end

Example
=======

given the ODE dy/dx = -2*y, x0=0, y0=1, xf=1, h=0.1, tolerance=1e-6. Use matrix M1 to store the rows of (x, y) values, skipping every 9 values.

The command is:

M1:=TRODE(0,1,1,0.1,0.000001,M1,9)

The value of y at xf (te last row of matrix M1) is calculated as 0.135335193015. The exact value is obtained using the command:

fxExact(1)

Which yields the value 0.135335283237.


Program Listing
===============

Code:

EXPORT fx(x,y)
BEGIN
  RETURN -2*y;
END;

// optional
EXPORT fxExact(x)
BEGIN
  RETURN EXP(-2*x);
END;

// function to solve in each ODE step
F(x0,y0,x1,y1,h)
BEGIN
  RETURN y1 - y0 - h/2*(fx(x0,y0) + fx(x1,y1));
END;

EXPORT TRODE(x0,y0,xf,h,toler,resMat,skipRows)
BEGIN
  LOCAL hx, x1, y1, F0, Fp, Fm;
  LOCAL diff, i, max, row, counter;
  counter := skipRows;
  resMat:=MAKEMAT(0,1,1);
  row := 1;
  resMat(row,1):=x0;
  resMat(row,2):=y0; 
  WHILE x0 < xf DO
    x1 := x0 + h;
    y1 := y0 + h*fx(x0,y0); // initial guess for y1
    diff := 2*toler;
    i := 1;
    max := 100;
    // loop to solve for y1 using Newton's method
    WHILE ABS(diff) >= toler AND i <= max DO
      i := i + 1;
      hx := 0.01*(1 + ABS(y1));
      F0 := F(x0,y0,x1,y1,h);
      Fp := F(x0,y0,x1,y1+hx,h);
      Fm := F(x0,y0,x1,y1-hx,h);
      diff := 2*hx*F0/(Fp - Fm);
      y1 := y1 - diff;
    END;
    y0 := y1;
    x0 := x1;
    IF counter < 1 THEN
      row := row + 1;
      resMat(row,1):=x0;
      resMat(row,2):=y0;
      counter:=skipRows;
    ELSE
      counter := counter - 1;
    END;
  END;
  IF resMat(row,1) < xf THEN
    row := row + 1;
    resMat(row, 1):= xf;
    resMat(row, 2):= y1;
  END;
  RETURN resMat;
END;