!********************************************************* !* Calculate a limited development of a real function * !* f(x) at point x0 with step h (to order 5). * !* ----------------------------------------------------- * !* SAMPLE RUN: * !* * !* Limited development of f(x) at x0: * !* * !* Function to develop: 1/(1-x) * !* * !* At point x0= 0 * !* Value of step h= 0.1 * !* * !* a1 = 1.0000416250 * !* a2 = 1.0000416250 * !* a3 = 1.0011207928 * !* a4 = 1.0011207928 * !* a5 = 1.0136468470 * !* * !* (The correct answer is a1=a2=a3=a4=a5=1) * !* * !* Function to develop: Ln(x) * !* * !* At point x0= 1 * !* Value of step h= 0.1 * !* * !* a1 = 1.0000057552 * !* a2 = -0.5000050520 * !* a3 = 0.3334508142 * !* a4 = -0.2501062334 * !* a5 = 0.2011301593 * !* * !* (the correct answer is a1=1, a2=-0.5, a3=1/3, * !* a4=-1/4 and a5=1/5) * !* ----------------------------------------------------- * !* Ref.: "Mathematiques en Turbo Pascal By Alain * !* Reverchon and Marc Ducamp, Editions Eyrolles, * !* Paris, 1991" [BIBLI 03]. * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !********************************************************* !Note: as you can see this direct method is not suitable to ! calculate a limited development in an accurate way. ! The next programs (ltddev1 to ltddev3) will allow us ! to calculate limited developments at x=0 up to a high ! order with a good accuracy for functions f*g, f/g or ! gof, knowing the limited developments of f(x) and g(x) ! at x=0. !-----------------------------------------------------------} Program LTDDEV real*8 d,h,r,s,x0; integer n print *,' ' print *,' Limited development of a real function f(x):' print *,' ' print *,' Function to develop: 1/(1-x)' print *,' ' write(*,10,advance='no'); read *, x0 write(*,20,advance='no'); read *, h print *,' ' s=1.d0 do n=1, 5 s = s * n r = ar_nderiv(x0,n,h) / s write(*,50) n, r end do print *,' ' 10 format(' At point x0= ') 20 format(' Value of step h= ') 50 format(' a',I1,' = ',F13.10) END !sample function f(x)=1/(1-x) real*8 Function f(x) real*8 x if (dabs(x-1) >1.d-16) then f = 1.d0/(1.d0-x) else f = 1.d16 end if end !sample function f(x)=Ln(x) !real*8 Function f(x) !real*8 x ! if (x >1.d-16) then ! f = LOG(x) ! else ! f = -1.d16 !end } !returns the value of nth derivative of function !f(x) at point x0 with increment step h (n=1 to 5). !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 x,h integer table(1:5,0:6) real*8 f,fa,fb,fc,t DATA table /3,3,4,4,5,45,270,-488,-1952,1938, & -9,-27,338,676,-1872,1,2,-72,-96,783, & 0,0,7,7,-152,0,0,0,0,13, & 60,180,240,240,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 End !end of file ltddev.f90