! ! Numerical Mathematics and Computing, Fifth Edition ! Ward Cheney & David Kincaid ! Brooks/Cole Publ. Co. ! Copyright (c) 2003. All rights reserved. ! For educational use with the Cheney-Kincaid textbook. ! Absolutely no warranty implied or expressed. ! ! Section 3.1 ! ! File: coef.f90 ! ! Newton interpolation polynomial for sin(x) at equidistant points (coef,eval) ! 10 equidistant points (9 intervals) ! evaluate sin(x)-p(x) at 37 points ! program main integer, parameter :: n=9 real, dimension (0:n) :: x, y, a real :: eval integer :: i interface function eval(n,x,a,t) integer :: n real, dimension(0:n) :: x, a real :: t end function eval subroutine coef(n,x,y,a) integer, intent(in) :: n real, dimension(0:n), intent(in) :: x, y real, dimension(0:n), intent(out) :: a end subroutine coef end interface h = 1.6875/real(n) do i=0,n x(i) = real(i)*h y(i) = sin(x(i)) end do call coef(n,x,y,a) print *, "Print the array a" print *, a print *, " Evaluation of sin(x)-p(x): n t d" do i = 0,4*n t = real(i)*h*0.25 d = sin(t) - eval(n,x,a,t) print *, i,t,d end do end program main subroutine coef(n,x,y,a) integer, intent(in) :: n integer :: i, j real, dimension(0:n), intent(in) :: x, y real, dimension(0:n), intent(out) :: a do i = 0,n a(i) = y(i) end do do j = 1,n do i = n,j,-1 a(i) = (a(i) - a(i-1))/(x(i) - x(i-j)) end do end do end subroutine coef function eval(n,x,a,t) integer :: i, n real, dimension(0:n) :: x, a real :: t, temp temp = a(n) do i=n-1,0,-1 temp = temp*(t - x(i)) + a(i) end do eval = temp end function eval