program main !======================================================================= ! Chapter 9, code 1: adaptive RKF45 for one first-order ODE x' = f(t,x). ! Based on Ode1_Adaptive.m by Alexander Godunov (March 2022). ! ! The driver supplies initial data and opens the data file. The RKF45 ! subroutine contains the stages and ALL adaptive step decisions. ! It prints and saves the initial point and every accepted point. ! ! Example: x' = 100/(1+10000*t**2), x(-1) = atan(-100). ! Exact solution: x(t) = atan(100*t). The rapid change near t = 0 ! makes the reduction and subsequent growth of the step size visible. ! ! To solve another equation, change the initial data and function f. ! Set check_exact = .false. unless you also change the exact formula. !======================================================================= implicit none double precision t0, x0, t_end, h_first, tolerance, xf, exact double precision f integer accepted, rejected integer, parameter :: data_unit = 20 logical check_exact external f !-------------------- Initial data: edit here ---------------------------- t0 = -1.0d0 x0 = atan(100.0d0*t0) t_end = 1.0d0 h_first = 0.1d0 tolerance = 1.0d-6 ! Absolute local error tolerance check_exact = .true. !------------------------------------------------------------------------ open(unit=data_unit, file='ch09_rkf45_scalar.dat', & status='replace', action='write') write(*,'(a)') 'Adaptive Runge-Kutta-Fehlberg RKF45' write(*,'(a,es12.4)') 'Tolerance: ', tolerance write(*,'(a)') ' t numerical x h' write(data_unit,'(a)') '# t x h' call rkf45(f, t0, x0, t_end, h_first, tolerance, data_unit, & xf, accepted, rejected) close(data_unit) write(*,'(a,i8)') 'Accepted steps: ', accepted write(*,'(a,i8)') 'Rejected trials:', rejected write(*,'(a,f18.12)') 'Final x: ', xf if (check_exact) then exact = atan(100.0d0*t_end) ! Exact formula for the DEFAULT example write(*,'(a,f18.12)') 'Exact final x: ', exact write(*,'(a,es12.4)') 'Abs error: ', abs(xf-exact) end if write(*,'(a)') 'Data file: ch09_rkf45_scalar.dat' end program main function f(t,x) !======================================================================= ! Supply x' = f(t,x) here. This example depends only on t, but the ! same subroutine also handles equations that depend on both t and x. !======================================================================= implicit none double precision f, t, x f = 100.0d0/(1.0d0+10000.0d0*t*t) end function f subroutine rkf45(f,t0,x0,t_end,h_first,tolerance,data_unit, & xf,accepted,rejected) !======================================================================= ! Integrate one first-order ODE with adaptive Fehlberg RKF45. ! Input: ! f derivative function f(t,x), supplied by the user ! t0, x0 initial time and value ! t_end requested final time (greater than t0) ! h_first first trial step size ! tolerance positive absolute local error tolerance ! data_unit unit of the data file already opened by the driver ! Output: ! xf solution at t_end ! accepted number of accepted steps (excluding the initial point) ! rejected number of rejected trials ! ! Six stages give fourth- and fifth-order estimates. Advance with ! the FIFTH-order value. Their difference estimates local error; ! tolerance is not a bound on total accumulated error. ! A rejected trial leaves both t and x unchanged. !======================================================================= implicit none double precision f, t0, x0, t_end, h_first, tolerance, xf double precision t, x, tf, h, trial, err, factor double precision k1, k2, k3, k4, k5, k6 integer data_unit, accepted, rejected, attempts integer, parameter :: max_attempts = 100000 double precision, parameter :: h_min = 1.0d-12 external f if (t_end <= t0 .or. h_first <= 0.0d0 .or. tolerance <= 0.0d0) & stop 'Use t_end > t0, h_first > 0, and tolerance > 0.' t = t0 x = x0 h = h_first accepted = 0 rejected = 0 attempts = 0 write(*,100) t, x, 0.0d0 write(data_unit,101) t, x, 0.0d0 do while (t < t_end) ! Ignore only a machine-rounding gap at the endpoint. if (t_end-t <= 64.0d0*epsilon(t)*max(1.0d0,abs(t_end))) exit if (attempts >= max_attempts) stop 'Maximum number of trials exceeded.' ! Shorten the trial if needed to finish at the requested time. h = min(h,t_end-t) tf = t+h if (h < h_min .or. tf <= t) stop 'Step too small to advance time.' attempts = attempts+1 ! Fehlberg stages: each k includes the factor h. k1 = h*f(t,x) k2 = h*f(t+h/4.0d0,x+k1/4.0d0) k3 = h*f(t+3.0d0*h/8.0d0,x+3.0d0*k1/32.0d0+9.0d0*k2/32.0d0) k4 = h*f(t+12.0d0*h/13.0d0,x+1932.0d0*k1/2197.0d0 & -7200.0d0*k2/2197.0d0+7296.0d0*k3/2197.0d0) k5 = h*f(t+h,x+439.0d0*k1/216.0d0-8.0d0*k2 & +3680.0d0*k3/513.0d0-845.0d0*k4/4104.0d0) k6 = h*f(t+h/2.0d0,x-8.0d0*k1/27.0d0+2.0d0*k2 & -3544.0d0*k3/2565.0d0+1859.0d0*k4/4104.0d0 & -11.0d0*k5/40.0d0) ! Fifth-order candidate and absolute fourth/fifth-order difference. trial = x+16.0d0*k1/135.0d0+6656.0d0*k3/12825.0d0 & +28561.0d0*k4/56430.0d0-9.0d0*k5/50.0d0+2.0d0*k6/55.0d0 err = abs(k1/360.0d0-128.0d0*k3/4275.0d0 & -2197.0d0*k4/75240.0d0+k5/50.0d0+2.0d0*k6/55.0d0) if (err <= tolerance) then if (h == t_end-t) tf = t_end t = tf x = trial accepted = accepted+1 ! h is the step used to reach this accepted point. write(*,100) t, x, h write(data_unit,101) t, x, h else rejected = rejected+1 end if ! The embedded difference scales approximately as h**5. ! Use a safety factor of 0.9 and limit growth/shrinkage as in ! the original MATLAB code. Zero error allows doubling of h. if (err == 0.0d0) then factor = 2.0d0 else factor = 0.9d0*(tolerance/err)**0.2d0 end if if (err <= tolerance) then h = h*min(2.0d0,factor) else h = h*max(0.2d0,factor) end if end do xf = x 100 format(f10.6,1x,f18.12,1x,es14.6) 101 format(3(1x,es24.16)) end subroutine rkf45