!
! =====> Program - P123.F90
!
PROGRAM DifDmo
! Program to demonstrate numerical solutions
! to ordinary, first order differential equations.
INTEGER I, N
REAL X, Xf, h, &
Y0, Y_Eul, Y_RK4
PRINT *, 'This is Program >> P123 = Differential Equations'
!
! Tell program where data for READ * is coming from
OPEN(UNIT=5, FILE='P123.DAT') ! UNIT=5 is the default input
!
PRINT *, 'First order ODEs F = -y*y for example'
PRINT *, 'Enter final (high) x value: '
READ *, Xf
Print *, Xf
PRINT *, 'Number of steps to take: '
READ *, N
Print *, N
! Step size, initial X, Y (given)
h = Xf/N
X = 0.0
Y0 = 1
PRINT *, 'Step size: ', h
PRINT 101, 'Step', 'x value', 'Euler', 'RK-4'
! Compare Euler and Runge-Kutta 4'th order
! at each step.
Y_Eul = Y0
Y_RK4 = Y0
L1: DO I = 1, N
Y_Eul = Euler(X, Y_Eul, h)
Y_RK4 = RngKt4(X, Y_RK4, h)
X = X + h
PRINT 103, I, X, Y_Eul, Y_RK4
END DO L1
STOP
101 FORMAT(//1X, A5, 6A12)
103 FORMAT(1X, I5, 3F12.8)
END
FUNCTION F(x, y)
! The function to be integrated.
! I.e. y' = dy/dx = F(x, y)
F = x ! Just to keep FTN90 quiet
F = -y*y
RETURN
END
! ALGORITHMS
FUNCTION Euler(x, y, h)
! Uses F(x, y) to advance the solution
! from x to x+h.
REAL x, y, h
Euler = y + h*F(x, y)
RETURN
END
FUNCTION RngKt4(x, y, h)
! Uses F(x, y) to advance the solution
! from x to x+h.
REAL x, y, h, &
h_hf, &
x_hf, &
y_hf_1, &
y_hf_2, &
y_1
h_hf = h/2
x_hf = x + h_hf
! Form various estimates of the solution
! at x+h/2 and x+h, then combine them
! using the formula below.
y_hf_1 = Euler(x, y, h_hf)
y_hf_2 = Euler(x_hf, y_hf_1, h_hf)
y_1 = Euler(x_hf, y_hf_2, h)
RngKt4 = y + h/3 * &
((f(x, y) + f(x+h, y_1))/2 + &
f(x_hf, y_hf_1) + f(x_hf, y_hf_2))
RETURN
END
DATA:
2.0
10
OUTPUT:
Program entered
This is Program >> P123 = Differential Equations
First order ODEs F = -y*y for example
Enter final (high) x value:
2.0000000
Number of steps to take:
10
Step size: 0.2000000
Step x value Euler RK-4
1 0.20000000 0.80000001 0.85231536
2 0.40000000 0.67199999 0.74049305
3 0.60000002 0.58168321 0.65333771
4 0.80000001 0.51401215 0.58374416
5 1.00000000 0.46117046 0.52703171
6 1.20000005 0.41863483 0.48001244
7 1.40000010 0.38358381 0.44045072
8 1.60000014 0.35415649 0.40673780
9 1.80000019 0.32907113 0.37768960
10 2.00000024 0.30741357 0.35241693
Fortran-90 STOP
Page builder: Charles Boivin
Last modified: 11/07/95