!********************************************************* !* Estimate the nth derivative of a real function f(x) * !* at point x with step h (n from 1 to 5). * !* ----------------------------------------------------- * !* SAMPLE RUN: * !* * !* Nth derivative of a real function f(x): * !* * !* Function to derivate: x^3 * !* * !* At point x0= 10 * !* Order of derivative (1 to 5): 1 * !* Value of step h= 0.1 * !* * !* f'(x0)=300.000000 * !* * !* At point x0= 10 * !* Order of derivative (1 to 5): 2 * !* Value of step h= 0.1 * !* * !* f"(x0)= 60.000000 * !* * !* ----------------------------------------------------- * !* Ref.: "Analyse numérique en C By Alain Reverchon and * !* Marc Ducamp, Armand Colin Editeur, Paris, * !* 1990" [BIBLI 14]. * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !********************************************************* Program Nderiv real*8 d,h,x integer n print *,' ' print *,' Nth derivative of a real function f(x):' print *,' ' print *,' Function to derivate: x^3' print *,' ' write(*,10,advance='no'); read *, x write(*,20,advance='no'); read *, n write(*,30,advance='no'); read *, h print *,' ' d = ar_nderiv(x,n,h) Select Case (n) case (1) write(*,50) d case (2) write(*,51) d case (3) write(*,52) d case (4) write(*,53) d case (5) write(*,54) d End Select if (n<1.or.n>5) & print *,' n must be between 1 and 5 !' print *,' ' 10 format(' At point x0= ') 20 format(' Order of derivative (1 to 5): ') 30 format(' Value of step h= ') 50 format(' f''(x0)=',F10.6) 51 format(' f"(x0)=',F10.6) 52 format(' f"''(x0)=',F10.6) 53 format(' f""(x0)=',F10.6) 54 format(' f""''(x0)=',F10.6) stop END !sample function f(x)=x^3 real*8 Function f(x) real*8 x f = x*x*x return end ! returns the value of nth derivative of function ! f(x) at point x0 with increment step h ! Note that for polynomial functions h must not be ! too small. Here h=1 gives best results than h=0.01 real*8 Function ar_nderiv(x, n, h) real*8 f,fa,fb,fc,fd,h,t,x integer n, i integer table(1:5,0:6) table(1,0)=3; table(1,1)=45; table(1,2)=-9; table(1,3)=1 table(1,4)=0; table(1,5)=0; table(1,6)=60 table(2,0)=3; table(2,1)=270; table(2,2)=-27; table(2,3)=2 table(2,4)=0; table(2,5)=0; table(2,6)=180 table(3,0)=4; table(3,1)=-488; table(3,2)=338; table(3,3)=-72 table(3,4)=7; table(3,5)=0; table(3,6)=240 table(4,0)=4; table(4,1)=-1952; table(4,2)=676; table(4,3)=-96 table(4,4)=7; table(4,5)=0; table(4,6)=240 table(5,0)=5; table(5,1)=1938; table(5,2)=-1872; table(5,3)=783 table(5,4)=-152; table(5,5)=13; table(5,6)=288 if (n>5) then ar_nderiv=0.d0 return end if fa = f(x) t=0.d0 do i=1, table(n,0) fb = f(x-i*h) fc = f(x+i*h) if (MOD(n,2).ne.0) then t = t + table(n,i)*(fc-fb) else t = t + table(n,i)*((fc-fa) + (fb-fa)) end if end do t = t / table(n,6) do i=1, n t = t / h end do ar_nderiv= t return End !end of file nderiv.f90