Fortran Need help with fortran 90 on an exercice

[fluidistic]

It's about Numerical Analysis and especially Newton's form of interpolation. That is, the program must ask us to input the nodes where the function (cos(x) in this case) is equal to it's interpolating polynomial.
From that P(x)=a_0+a_1(x-x_0)+...+a_n(x-x_o)...(x-x_{n-1}).
And then it must output the value of the interpolating polynomial evaluated m times between x_0 and x_n. m is a value written in the program (that you can change), of about 100 for example. It means the output of the program will give 100 values of P(x) in the interval x_0, x_n. m divided the interval x_0, x_n properly (all the subintervals have same length).
So I worked out my program... But when I graph the output via Gnuplot... It doesn't look like the cos(x) function at all! I would like to see where is my error(s).
If someone could do this, I would be very glad. If you have some doubts, just ask.
Here comes the program... :

Program Interp
implicit none

Integer, Parameter :: n=5
Real,Dimension(0:n-1) :: x,y
Integer :: K,J,m
Real :: z,p,dx,s

m=100

do k=0,n-1
Write(*,*)'Enter x sub',K
Read(*,*)x(K)
y(K)=f(x(K))
end do

dx=(x(n-1)-x(0))/m
z=x(0)
DO WHILE (z<=x(n-1))
CALL NEWTON(z,p,x,y,s)
WRITE(*,*)z,s
z=z+dx
end doContains
Subroutine Newton(z,p,x,y,s)
implicit none
Real, Dimension(0:n-1),Intent(in) ::x,y
Real, Intent(out) :: S
Real, Intent(in) :: z
Real, Dimension(0:n-1) :: a
Integer :: K,J,L
Real :: D,P

a(0)=y(0)
a(1)=(y(1)-a(0))/(x(1)-x(0))
Do K=2,n-1,1
D=1
Do J=0,K-1
D=D*(x(k)-x(J))
end do
P=A(k-1)
Do L=K-1,0,-1
P=A(L)+P*(x(K)-x(L))
end do
A(K)=(Y(K)-P)/D
end do
S=A(n-1)
DO K=n-1,0,-1
S=A(K)+S*(z-x(K))
end do
end subroutine
Real Function f(x)
Real, intent(in) :: x
f=cos(x)
end function

end program

 

[alphysicist]

Hi fluidistic,

I believe I see two errors in your program. You calculate a_2 \to a_n using:

D=1
Do J=0,K-1
D=D*(x(k)-x(J))
end do
P=A(k-1)
Do L=K-1,0,-1 <-----------------------error?
P=A(L)+P*(x(K)-x(L))
end do
A(K)=(Y(K)-P)/D

From the polynomial we need, for example,

a_2 = \frac{y_2 -a_1(x_2 - x_0) - a_0}{(x_2-x_0)(x_2-x_1)}

To get this I think you need to first line of the do loop to be

do L=k-2,0,-1

or else you'll get two terms with a_{k-1} in your formula for a_k.
---------------------------------------------

I think the same error is present when you construct the polynomial S:

S=A(n-1)
DO K=n-1,0,-1 <-------------error?
S=A(K)+S*(z-x(K))
end do

With the do loop beginning at n-1 as it is here, the polynomial will have two terms with a_{n-1}. I think you want to have it be

DO K=n-2,0,-1


-------------------------

I ran the program using x values of {0,1,2,3,4} as well as some others and at least visually seemed to fit fine (if the spacing was not too great, of course).

 

[fluidistic]

Thanks a lot

Thank you very very much! It does work as it should now...
I'm learning how to program with Fortran and I find very hard some exercises. Anyway it wasn't a homework, it was a preparation for an exam I'll get on tuesday, so now I can study the program better. Don't know how to thank you!

 

[alphysicist]

I'm glad I could help! In my work for the past couple of weeks I've been caculating derivatives of tabulated functions and so these interpolation calculations have been on my mind a lot.

 

[ice109]

the Newton form of the interpolating polynomial is more difficult to code than the lagrange. i would suggest using the lagrange if you can.

 

[fluidistic]

Hello ice109,
In my anterior exercise I had to do it with Lagrange's interpolation and it worked. For me both are difficult!