! ! =====> 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