02-20-2024, 10:37 PM
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:
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
===============
================================================
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
endExample
=======
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;