PROGRAM Simple_Differential ! ! This program solves a simple differential equation given by ! ! y'' + y' = 0 ! ! with initial conditions such that ! ! y = 1 for x=0 & y'=-1 for x=0. ! ! The solution is given by y(x) = exp(-x). ! ! The code generates a text file, "solution.txt", which contains ! ! x, y(x) numerically solved, the analytical solution for y(x) ! ! You can check whether the numerical solution agrees with the exact ! solution. ! ! The differntial equation is solved by a publicly available routine ! called "dverk.f". For more information on dverk.f, see ! http://www.netlib.org/ode/dverk.f ! ! E.Komatsu, March 31, 2008 ! IMPLICIT none DOUBLE PRECISION :: x,xinitial,y(2),c(24),work(2,9),tol=1d-8 INTEGER :: i,ind EXTERNAL fderiv open(1,file='solution.txt') ind=1 do i=0,500 x = i*0.01d0 if (i.eq.0) then xinitial = x y(1) = 1d0 ! set initial condition, such that y=1 for x=0 y(2) = -1d0 ! set initial condition, such that y'=-1 for x=0 else call dverk(2,fderiv,xinitial,y,x,tol,ind,c,2,work) xinitial = x endif write(1,'(3E20.8)')x,y(1),dexp(-x) enddo close(1) stop END PROGRAM Simple_Differential !----------------------------------------------------------------------------- SUBROUTINE fderiv(n,x,y,yprime) ! ! Define a differential equation given by ! ! d^2y/dx^2 + dy/dx = 0 ! ! whose solution is ! ! y(x) = a*exp(-x) + b ! ! where a and b are integration constants. ! ! In this subroutine, we need to define ! - x ! - y(1) = y(x) ! - y(2) = dy/dx ! IMPLICIT none INTEGER, intent(IN) :: n ! # of derivatives DOUBLE PRECISION :: x,y(n),yprime(n) yprime(1)= y(2) ! y(1) is y. So, "yprime(1)" is y' = dy/dx yprime(2)= -y(2) ! y(2) is dy/dx. So, "yprime(2)" is y'' = d^y/dx^2 return END SUBROUTINE fderiv