hermite.for

一个非常好的插值程序

	program main
	double precision  x(10),y(10),dy(10),x1,y1
	data x/0.10,0.15,0.30,0.45,0.55,0.60,0.70,0.85,
	1      0.90,1.0/
	data y/0.904837,0.860708,0.740818,0.637628,0.576950,
	1      0.548812,0.496585,0.427415,0.406570,0.367879/
	do 1 i=1,10
	  dy(i)=-y(i)
1	continue
	call hermite(x,y,dy,10,x1,y1)
	end

      subroutine hermite(x,y,dy,n,x1,y1)
	double precision x(n),y(n),dy(n),x1,y1,s,p,q

	y1=0.0		 !good
	print *,"please input x value!"
	read *,x1

	if((x1-x(1)).lt.-1e-6.or.x1.gt.x(n)) then
	  print *,"input value is not in the content!"
	else
	  if(n.lt.0) then
	  return
	  else if(n.eq.1) then
	  y1=y(1)
	  else if(n.eq.2) then
	  y1=(y(1)*(x1-x(2))-y(2)*(x1-x(1)))/(x(1)-x(2))
	  else
	      i=1
10	      if(x(i).lt.x1) then
	         i=i+1
	         if(i.le.n) goto 10
	      end if
	      k=i-4
	      if(k.lt.1) k=1
	      m=i+3
	      if(m.gt.n) m=n
	  end if
	  do 30 i=k,m
	    s=1.0			!good
	    do 20 j=k,m
	      if(j.ne.i) then
	        s=s*(x1-x(j))/(x(i)-x(j))
	      end if
20	    continue
          s=s**2
	    p=0.0
	    do 25 j=k,m
	       if(j.ne.i) then
	         p=p+1.0/(x(i)-x(j))
	       endif
25	    continue
          q=y(i)+(x1-x(i))*(dy(i)-2*y(i)*p)
	    y1=y1+s*q
30	  continue
	  print *,x1,y1
	end if
      end