newton.for

非常好的插值程序

      	program main
	double precision  x(10),y(10),x1,y1
	data x/0.10,0.15,0.25,0.40,0.50,0.57,0.70,0.85,
	1      0.93,1.0/
	data y/0.904837,0.860708,0.778801,0.670320,0.606531,
	1      0.565525,0.496585,0.427415,0.394554,0.367879/
	call newton(x,y,10,x1,y1)
	end

      subroutine newton(x,y,n,x1,y1)
	double precision x(n),y(n),x1,y1,b(n),z,u,s

	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
	      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 j=k+1,m
	    u=0.0      !good
	   do 70 i=k,j
            s=1.0     !good
	      do 20 l=k,j
		   	 if(l.ne.i) then	          
	           s=s*1.0/(x(i)-x(l))
	         endif
20          continue
            u=u+y(i)*s
70	   continue                    
          b(j)=u
30    continue
 
      y1=y(k)	          !good
	do 50 j=k+1,m
	  z=1.0	            !good
        do 40 i=k,j-1
	     z=z*(x1-x(i))
40      continue
       y1=y1+z*b(j)
50	continue
	print *,x1,y1
	end if
      end