Implicit Trapezoidal Integration Method for ODE
02-20-2024, 10:37 PM (This post was last modified: 02-21-2024 12:40 PM by Namir.)
Post: #1
 Namir Senior Member Posts: 1,030 Joined: Dec 2013
Implicit Trapezoidal Integration Method for ODE
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;
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)