c-----*----------------------------------------------------6---------7-- c This subroutine is to calculate the m-th numerical derivative of c unequally spaced data by use of an (n+1)-point formula. c Here the maximum of m is 8. c You can easily compile a subroutine for higher derivative than 8 through c modifying the subroutine root(n,m,x,a). c input: nn,n,m,x,f c nn: (nn+1) is the number of sampling data, from 0 to nn c n: (n+1)-point formula,from 0 to n, n< or =nn c m: m-th derivative, m is less than or equal to n. c x(0:nn): (nn+1) sampling data c f(0:nn): function values at (nn+1) points c output: df c df(0:nn): approximations of the m-th derivative at (nn+1) sampling points c work array: c(0:n,0:n) c c(0:n,0:n): stores (n+1) coefficients of the (n+1)-point formula, c In d(j,i), the index i=0,n denote (n+1) points xi (i=0,n), c and the index j=0,n represent (n+1) coefficients of m-th numerical c derivative at the conference point xi (i=0,n). c subroutine: uneqcoeff is to obtain the the coefficients d of unequally c spaced data. cccccc c Reference: c Li, J., 2005: General explicit difference formulas for numerical differentiation. c J. Comput. Appl. Math., 183(1), 29-52, doi: 10.1016/j.cam.2004.12.026 cccccc c By Dr. Jianping Li, April 2, 2003. subroutine Udiffer(nn,n,m,x,f,df) dimension x(0:nn),f(0:nn),df(0:nn) real y(0:n),c(0:n) !work array do i=0,nn df(i)=0. enddo c---- k=int(n/2.) do i=0,k-1 do j=0,n y(j)=x(j) enddo call uneqcoeff(n,m,i,y,c) do j=0,n if(j.ne.i)then d=(f(j)-f(i))/(x(j)-x(i)) df(i)=df(i)+d*c(j) endif enddo enddo do i=k,nn-k-1 ik=i-k do j=0,n y(j)=x(ik+j) enddo call uneqcoeff(n,m,k,y,c) do j=0,n jk=ik+j if(jk.ne.i)then d=(f(jk)-f(i))/(x(jk)-x(i)) df(i)=df(i)+d*c(j) endif enddo enddo nk=nn-n do i=nn-k,nn ii=n-(nn-i) do j=0,n y(j)=x(nk+j) enddo call uneqcoeff(n,m,ii,y,c) do j=0,n jk=nk+j if(jk.ne.i)then d=(f(jk)-f(i))/(x(jk)-x(i)) df(i)=df(i)+d*c(j) endif enddo enddo return end c-----*----------------------------------------------------6---------7-- c This subroutine is used to calculate coefficients of an (n+1)-point c formula for the m-th numerical derivative of unequally spaced data. c Here the maximum of m is 8. c One can easily calculate those coefficients for larger m c through inputting the common denominator array and c writing the corresponding code in the subroutine ROOT. c input: n,m,i,x c n: (n+1)-point formula c m: m-th derivative, m is less than or equal to n. c i: the given conference point c x(0:n): (n+1) numbers c output: c c c(0:n): coefficient array for the conference point xi c work array: x c work function: ifact c ifact(i): external function, factorial of i c By Dr. Jianping Li, April 2, 2003. subroutine uneqcoeff(n,m,i,x,c) dimension c(0:n),x(0:n) dimension y(n),z(n) !work array real ifact If(m.gt.n)then write(*,*)'m must be less than or equal to n' stop endif fm=ifact(m) do j=0,n c(j)=0. enddo do j=0,n if(j.ne.i)then kn=0 do k=0,n if(k.ne.i.and.k.ne.j)then kn=kn+1 y(kn)=x(k)-x(i) z(kn)=x(k)-x(j) endif enddo call root(kn,m-1,y,a) c1=(-1)**(m-1)*fm*a call rfact(kn,z,c2) c(j)=c1/c2 endif enddo return end c-----*----------------------------------------------------6---------7-- c Computing the factorial n! of a nutural number n. real function ifact(n) if(n.lt.0)then print*,'Illegal factorial',n stop else ifact=1.0 do i=2,n ifact=ifact*float(i) enddo endif return end c-----*----------------------------------------------------6---------7-- c This subroutine is to calculate the root coefficient of x1,...,xn c (m) c An =Am(x1,x2,...,xn). c input: n,m,x c n: number of x c m: superscript index of An(m) c x(n): n numbers c output: a c a: =An(m)=Am(x1,x2,...,xn) c By Dr. Jianping Li, April 2, 2003. subroutine root(n,m,x,a) dimension x(n) if(m.eq.0)then a=1. do i=1,n a=a*x(i) enddo return endif if(m.eq.1)then a=0. do j=1,n b=1. do i=1,n if(i.ne.j)then b=b*x(i) endif enddo a=a+b enddo return endif if(m.eq.2)then a=0. do k=1,n-1 do j=k+1,n b=1. do i=1,n if(i.ne.j.and.i.ne.k)then b=b*x(i) endif enddo a=a+b num=num+1 enddo enddo return endif if(m.eq.3)then a=0. do l=1,n-2 do k=l+1,n-1 do j=k+1,n b=1. do i=1,n if(i.ne.j.and.i.ne.k.and.i.ne.l)then b=b*x(i) endif enddo a=a+b enddo enddo enddo return endif if(m.eq.4)then a=0. do 40 k1= 1,n-3 do 40 k2=k1+1,n-2 do 40 k3=k2+1,n-1 do 40 k4=k3+1,n b=1. do i=1,n if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4)then b=b*x(i) endif enddo a=a+b 40 continue return endif if(m.eq.5)then a=0. do 50 k1= 1,n-4 do 50 k2=k1+1,n-3 do 50 k3=k2+1,n-2 do 50 k4=k3+1,n-1 do 50 k5=k4+1,n b=1. do i=1,n if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and. * i.ne.k5)then b=b*x(i) endif enddo a=a+b 50 continue return endif if(m.eq.6)then a=0. do 60 k1= 1,n-5 do 60 k2=k1+1,n-4 do 60 k3=k2+1,n-3 do 60 k4=k3+1,n-2 do 60 k5=k4+1,n-1 do 60 k6=k5+1,n b=1. do i=1,n if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and. * i.ne.k5.and.i.ne.k6)then b=b*x(i) endif enddo a=a+b 60 continue return endif if(m.eq.7)then a=0. do 70 k1= 1,n-6 do 70 k2=k1+1,n-5 do 70 k3=k2+1,n-4 do 70 k4=k3+1,n-3 do 70 k5=k4+1,n-2 do 70 k6=k5+1,n-1 do 70 k7=k6+1,n b=1. do i=1,n if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and. * i.ne.k5.and.i.ne.k6.and.i.ne.k7)then b=b*x(i) endif enddo a=a+b 70 continue return endif if(m.eq.8)then a=0. do 80 k1= 1,n-7 do 80 k2=k1+1,n-6 do 80 k3=k2+1,n-5 do 80 k4=k3+1,n-4 do 80 k5=k4+1,n-3 do 80 k6=k5+1,n-2 do 80 k7=k6+1,n-1 do 80 k8=k7+1,n b=1. do i=1,n if(i.ne.k1.and.i.ne.k2.and.i.ne.k3.and.i.ne.k4.and. * i.ne.k5.and.i.ne.k6.and.i.ne.k7.and.i.ne.k8)then b=b*x(i) endif enddo a=a+b 80 continue return endif return end c-----*----------------------------------------------------6---------7-- c Computing the product of n arbitrary numbers. c input: n, x(n) c output: p c p=x(1)*x(2)*...*x(n) subroutine rfact(n,x,p) dimension x(n) p=1. do i=1,n p=p*x(i) enddo return end