! ! =====> Program - P123.F90 ! PROGRAM DifDmo IMPLICIT NONE ! 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' ! ! 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(//' ', A5, 6A12) 103 FORMAT(' ', I5, 3F12.8) END PROGRAM DifDmo 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 FUNCTION F ! 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 Euler 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 FUNCTION RngKt4