Fortran, interpolating polynomial program doesn't compile

[fluidistic]

I want to interpolate a set of points (x(i),y(i)) for i=0, n; using Newton's polynomial of degree n.
I've Kincaid's book of numerical analysis. He gives an algorithm that doesn't make entirely sense to me.
He says for

i=0, n do 
di <- f(xi)
end
for j=1, n do 
etc.

Basically I'm stuck in Fortran for the di <- f(xi) part.
I'll show you my whole program. What I did doesn't make sense to gfortran (the compiler), but I do not understand why!

program div
implicit none
character(10000000) :: datos
real, dimension(:), allocatable :: x,y,d,p
integer::n ,i
real :: h, x_m
write(*,*)"Enter the degree of the polynomial n and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file to read the data from"
read(*,'(x)')datos
open(unit=1,file="datos", status="old")
allocate(x(0:n),y(0:n), d(1:n), p(1:n))
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)

call Newton(x_m,h,x,d,y)

contains
subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h,d
real:: prod, p
real, dimension(:), intent(inout) :: x,y

open(unit=120, file="polynomial.dat")
do while (x_m <= x(n))
  [b]do i=0,n
  d(i)=y(i)
  end do[/b]
            do j=1, n
            if (j>= i) then
            [b]d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))[/b]
            end if
            end do
prod=x_m-x(0)
                                            do i=1, n  
                                               do j=0, i-1
                                               prod=prod*x_m-x(j)
                                               end do
                                            p=d(i)*prod
                                            write(120,*)p
                                            end do    

x_m=x_m+h
end do
close(120)
end subroutine
end program

Here are the errors I'm getting:

Newton.f90:34.2:

  d(i)=y(i)
  1
Error: Unclassifiable statement at (1)
Newton.f90:38.12:

            d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
            1
Error: Unclassifiable statement at (1)
Newton.f90:24.27:

subroutine Newton(x_m,h,x,d,y)
                           1
Error: PROCEDURE attribute conflicts with INTENT attribute in 'd' at (1)
Newton.f90:41.11:

prod=x_m-x(0)
           1
Warning: Array reference at (1) is out of bounds (0 < 1) in dimension 1
Newton.f90:19.20:

call Newton(x_m,h,x,d,y)
                    1
Error: Rank mismatch in argument 'd' at (1) (0 and 1)

Some more details about the program: I have to write the name of the file I read the points from, say "kincaid.dat" and the degree of the polynomial. For example if I have 10 points to interpolate, I'd enter n=9.
I want my program to write in a file (polynomial.dat) the values for p(x_m). If I enter a small h, I'll get lots of values for p(x_m) so that if I graph p(x_m) in gnuplot I can have a good idea of the theoretical graph of p(x), the theoretical interpolating polynomial of the data.
Any help is greatly appreciated. I feel like I'm turning crazy passing hours in a row on this program and now I'm stuck.

 

[Mark44]

I want to interpolate a set of points (x(i),y(i)) for i=0, n; using Newton's polynomial of degree n.
I've Kincaid's book of numerical analysis. He gives an algorithm that doesn't make entirely sense to me.
He says for

i=0, n do 
di <- f(xi)
end
for j=1, n do 
etc.

Basically I'm stuck in Fortran for the di <- f(xi) part.
I'll show you my whole program. What I did doesn't make sense to gfortran (the compiler), but I do not understand why!

Thanks for using code tags! I also fiddled with your indentation to make your code more readable.

program div
implicit none
character(10000000) :: datos
real, dimension(:), allocatable :: x,y,d,p
integer::n ,i
real :: h, x_m
write(*,*)"Enter the degree of the polynomial n and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file to read the data from"
read(*,'(x)')datos
open(unit=1,file="datos", status="old")
allocate(x(0:n),y(0:n), d(1:n), p(1:n))
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)

call Newton(x_m,h,x,d,y)

contains
subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h,d
real:: prod, p
real, dimension(:), intent(inout) :: x,y

open(unit=120, file="polynomial.dat")
do while (x_m <= x(n))
  [b]do i=0,n
     d(i)=y(i)
  end do[/b]
  do j=1, n
     if (j>= i) then
    [b]d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))[/b]
     end if
  end do
  prod=x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*x_m-x(j)
    end do
    p=d(i)*prod
    write(120,*)p
  end do    

  x_m=x_m+h
end do
close(120)
end subroutine
end program

Here are the errors I'm getting:

Newton.f90:34.2:

  d(i)=y(i)
  1
Error: Unclassifiable statement at (1)
Newton.f90:38.12:

            d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
            1
Error: Unclassifiable statement at (1)
Newton.f90:24.27:

subroutine Newton(x_m,h,x,d,y)
                           1
Error: PROCEDURE attribute conflicts with INTENT attribute in 'd' at (1)
Newton.f90:41.11:

prod=x_m-x(0)
           1
Warning: Array reference at (1) is out of bounds (0 < 1) in dimension 1
Newton.f90:19.20:

call Newton(x_m,h,x,d,y)
                    1
Error: Rank mismatch in argument 'd' at (1) (0 and 1)

Some more details about the program: I have to write the name of the file I read the points from, say "kincaid.dat" and the degree of the polynomial. For example if I have 10 points to interpolate, I'd enter n=9.
I want my program to write in a file (polynomial.dat) the values for p(x_m). If I enter a small h, I'll get lots of values for p(x_m) so that if I graph p(x_m) in gnuplot I can have a good idea of the theoretical graph of p(x), the theoretical interpolating polynomial of the data.
Any help is greatly appreciated. I feel like I'm turning crazy passing hours in a row on this program and now I'm stuck.

I believe that your problem stems from using d and x as arrays in your main program, but as scalars in your subroutine. In the subroutine d and x are not arrays, so the compiler is complaining about using these variables with indexes, as in d(i) and x(i).

 

[fluidistic]

Thanks for using code tags! I also fiddled with your indentation to make your code more readable.I believe that your problem stems from using d and x as arrays in your main program, but as scalars in your subroutine. In the subroutine d and x are not arrays, so the compiler is complaining about using these variables with indexes, as in d(i) and x(i).

Thanks for your help.
I tried to define new variables (scalars) in function of x(i), y(i) in the subroutine as to get rid of the arrays, but without success.
I'm really confused about how to solve the problem. Basically I must transform the arrays x and d into scalars, in the subroutine?

Edit: I tried to modify only a small part of the declaration of variables in the subroutine.
Now it looks like this:

subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d

The program compiles! But... I do get an error.
Here it is:

Enter the degree of the polynomial n and a step-size h
3, 0.001
 Enter the name of the file to read the data from
kincaid.dat
At line 10 of file hope2.f90 (unit = 5, file = 'stdin')
Fortran runtime error: Insufficient data descriptors in format after reversion
fluidistic@fluidistic-desktop:~/numerical2011/guide5$ ncaid.dat
/usr/bin/python: can't find '__main__' module in '/usr/share/command-not-found'

My line 10 is

read(*,'(x)')datos

I'm guessing the '(x)' part is a problem. I must say I randomly wrote '(x)'. My professor gave us a -working- program with such an example.
It starts like this:

character(50)::matrix, vector
  real, dimension(:,:), allocatable::A
  real, dimension(:), allocatable::b,x,order
  integer::n,i,j

  write(*,*)'Enter n'
  read(*,*)n

  write(*,*)'Enter the name of the file of the matrix A'
  read(*,'(A)')matrix
  open(10,file=matrix,status='old')

  write(*,*)'Do the same for the b vector'
  read(*,'(A)')vector
  open(15,file=vector,status='old')

  allocate(A(n,n),b(n),x(n),order(n))

But I really do not understand how the read(*,'whatisthis?') works.

Edit 2: I searched on the internet and I've learned that '(A)' has nothing to do with the A matrix. Instead, it's a format or something like that. So I replaced my '(x)' by '(A)'. I also erased some "".
Now my program looks like this:

program div
implicit none
character(100000) :: datos
real, dimension(:), allocatable :: x,y,d,p
integer::n ,i
real :: h, x_m
write(*,*)"Enter the degree of the polynomial n and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(unit=1,file=datos, status="old")
allocate(x(0:n),y(0:n), d(1:n), p(1:n))
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)

call Newton(x_m,h,x,d,y)

contains
subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=120, file="polynomial.dat")
do while (x_m <= x(n))
  do i=0,n
     d(i)=y(i)
  end do
  do j=1, n
     if (j>= i) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
  prod=x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*x_m-x(j)
    end do
    p=d(i)*prod
    write(120,*)p
  end do    

  x_m=x_m+h
end do
close(120)
end subroutine
end program

But I'm having an error of the form *** glibc detected *** ./a.out: free(): invalid next size (fast): 0x08f74018 ***.

I'm guessing it has to do with the '(A)' format?
I'm trying to read from a file with 2 columns of n rows. What should I use instead of '(A)'?

Edit3: I could compile+execute without error, though I had to modificate badly my program. I get only 1 number as output instead of more than 100.
I'll work on that...

 

[fluidistic]

It seems like I have at least 2 errors. I still get the memory error and the output (hundreds of lines) is wrong. The program can still run though. Could someone tell me what I'm doing wrong? I'm guessing I do not close a file when I should...
Here is my program:

program div
implicit none
character(100000) :: datos
real, dimension(:), allocatable :: x,y,d
integer::n ,i
real :: h, x_m
write(*,*)"Enter the degree of the polynomial n and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(1,file=datos, status='old')
allocate(x(0:n),y(0:n), d(1:n))
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)
x_m=x(0)
call Newton(x_m,h,x,d,y)

contains
subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=120, file="polynomial.dat", status='old')
do while (x_m <= x(n))
  do i=0,n
     d(i)=y(i)
  end do
  do j=1, n
     if (j>= i) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
  prod=x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*x_m-x(j)
    end do
    p=d(i)*prod
    write(120,*)x_m, p

  end do    

  x_m=x_m+h
end do
close(120)
end subroutine
end program

The name of the file I must input is "kincaid.dat" and the data inside the file is:

1 -3
3 1
5 2
6 4

Since it's 4 points, I must enter n=3 as input. I use a polynomial of degree 3 to interpolate those points.
My output file "polynomial.dat" is: (I've had to cut it to enter here. x_m does go from 1 to 6 as it should)

 1.0000000      -1.0000000    
   1.0000000      -10.000000    
   1.0000000      -56.000000    
   1.0100000     -0.98989999    
   1.0100000      -10.039594    
   1.0100000      -56.888004    
   1.0200000     -0.97960001    
   1.0200000      -10.078352    
   1.0200000      -57.792053    
   1.0300000     -0.96910006    
   1.0300000      -10.116237    
   1.0300000      -58.712166    
   1.0400000     -0.95840007    
   1.0400000      -10.153211    
   1.0400000      -59.648361    
   1.0500000     -0.94750005    
   1.0500000      -10.189238    
   1.0500000      -60.600628    
   1.0599999     -0.93640006    
   1.0599999      -10.224277    
   1.0599999      -61.568951    
   1.0699999     -0.92510009    
   1.0699999      -10.258294    
   1.0699999      -62.553299    
   1.0799999     -0.91360009    
   1.0799999      -10.291245    
   1.0799999      -63.553600    
   1.0899999     -0.90190011    
   1.0899999      -10.323094    
   1.0899999      -64.569801    
   1.0999999     -0.89000010    
   1.0999999      -10.353800    
   1.0999999      -65.601807    
   1.1099999     -0.87790012    
   1.1099999      -10.383321    
   1.1099999      -66.649490    
   1.1199999     -0.86560017    
   1.1199999      -10.411617    
   1.1199999      -67.712730    
   1.1299999     -0.85310018    
   1.1299999      -10.438646    
   1.1299999      -68.791367    
   1.1399999     -0.84040016    
   1.1399999      -10.464367    
   1.1399999      -69.885223    
   1.1499999     -0.82750016    
   1.1499999      -10.488737    
   1.1499999      -70.994102    
   1.1599998     -0.81440020    
   1.1599998      -10.511713    
   1.1599998      -72.117767    
   1.1699998     -0.80110019    
   1.1699998      -10.533251    
   1.1699998      -73.255966    
   1.1799998     -0.78760022    
   1.1799998      -10.553308    
   1.1799998      -74.408432    
   1.1899998     -0.77390027    
   1.1899998      -10.571839    
   1.1899998      -75.574837    
   1.1999998     -0.76000029    
   1.1999998      -10.588799    
   1.1999998      -76.754868    
   1.2099998     -0.74590027    
   1.2099998      -10.604144    
   1.2099998      -77.948151    
   1.2199998     -0.73160028    
   1.2199998      -10.617826    
   1.2199998      -79.154289    
   1.2299998     -0.71710032    
   1.2299998      -10.629801    
   1.2299998      -80.372864    
   1.2399998     -0.70240033    
   1.2399998      -10.640020    
   1.2399998      -81.603409    
   1.2499998     -0.68750036    
   1.2499998      -10.648438    
   1.2499998      -82.845428    
   1.2599998     -0.67240036    
   1.2599998      -10.655005    
   1.2599998      -84.098404    
   1.2699997     -0.65710038    
   1.2699997      -10.659673    
   1.2699997      -85.361755    
   1.2799997     -0.64160043    
   1.2799997      -10.662395    
   1.2799997      -86.634895    
   1.2899997     -0.62590045    
   1.2899997      -10.663120    
   1.2899997      -87.917175    
   1.2999997     -0.61000043    
   1.2999997      -10.661799    
   1.2999997      -89.207909    
   1.3099997     -0.59390050    
   1.3099997      -10.658384    
   1.3099997      -90.506401    
   1.3199997     -0.57760048    
   1.3199997      -10.652821    
   1.3199997      -91.811852    
   1.3299997     -0.56110054    
   1.3299997      -10.645060    
   1.3299997      -93.123482    
   1.3399997     -0.54440057    
   1.3399997      -10.635050    
   1.3399997      -94.440430    
   1.3499997     -0.52750057    
   1.3499997      -10.622738    
   1.3499997      -95.761787    
   1.3599997     -0.51040059    
   1.3599997      -10.608072    
   1.3599997      -97.086632    
   1.3699996     -0.49310061    
   1.3699996      -10.591000    
   1.3699996      -98.413948    
   1.3799996     -0.47560063    
   1.3799996      -10.571465    
   1.3799996      -99.742691    
   1.3899996     -0.45790067    
   1.3899996      -10.549418    
   1.3899996      -101.07178    
   1.3999996     -0.44000068    
   1.3999996      -10.524801    
   1.3999996      -102.40005    
   1.4099996     -0.42190072    
   1.4099996      -10.497560    
   1.4099996      -103.72630    
   1.4199996     -0.40360072    
   1.4199996      -10.467639    
   1.4199996      -105.04927    
   1.4299996     -0.38510075    
   1.4299996      -10.434983    
   1.4299996      -106.36764    
   1.4399996     -0.36640078    
   1.4399996      -10.399535    
   1.4399996      -107.68003    
   1.4499996     -0.34750080    
   1.4499996      -10.361239    
   1.4499996      -108.98500    
   1.4599996     -0.32840085    
   1.4599996      -10.320037    
   1.4599996      -110.28105    
   1.4699996     -0.30910087    
   1.4699996      -10.275870    
   1.4699996      -111.56660    
   1.4799995     -0.28960091    
   1.4799995      -10.228683    
   1.4799995      -112.84005    
   1.4899995     -0.26990092    
   1.4899995      -10.178412    
   1.4899995      -114.09967    
   1.4999995     -0.25000095    
   1.4999995      -10.125003    
   1.4999995      -115.34370    
   1.5099995     -0.22990099    
   1.5099995      -10.068393    
   1.5099995      -116.57029    
   1.5199995     -0.20960101    
   1.5199995      -10.008523    
   1.5199995      -117.77753    
   1.5299995     -0.18910104    
   1.5299995      -9.9453316    
   1.5299995      -118.96346    
   1.5399995     -0.16840108    
   1.5399995      -9.8787584    
   1.5399995      -120.12599    
   1.5499995     -0.14750110    
   1.5499995      -9.8087416    
   1.5499995      -121.26297    
   1.5599995     -0.12640113    
   1.5599995      -9.7352180    
   1.5599995      -122.37219    
   1.5699995     -0.10510116    
   1.5699995      -9.6581259    
   1.5699995      -123.45133    
   1.5799994     -8.36011916E-02
   1.5799994      -9.5774031    
   1.5799994      -124.49804    
   1.5899994     -6.19012266E-02
   1.5899994      -9.4929838    
   1.5899994      -125.50980    
   1.5999994     -4.00012583E-02
   1.5999994      -9.4048052    
   1.5999994      -126.48407    
   1.6099994     -1.79012921E-02
   1.6099994      -9.3128023    
   1.6099994      -127.41818    
   1.6199994      4.39867564E-03
   1.6199994      -9.2169113    
   1.6199994      -128.30942    
   1.6299994      2.68986430E-02
   1.6299994      -9.1170645    
   1.6299994      -129.15489    
   1.6399994      4.95986082E-02
   1.6399994      -9.0131979    
   1.6399994      -129.95172    
   1.6499994      7.24985749E-02
   1.6499994      -8.9052448    
   1.6499994      -130.69684    
   1.6599994      9.55985412E-02
   1.6599994      -8.7931366    
   1.6599994      -131.38710    
   1.6699994      0.11889850    
   1.6699994      -8.6768074    
   1.6699994      -132.01930    
   1.6799994      0.14239848    
   1.6799994      -8.5561886    
   1.6799994      -132.59009    
   1.6899993      0.16609843    
   1.6899993      -8.4312124    
   1.6899993      -133.09599    
   1.6999993      0.18999840    
   1.6999993      -8.3018084    
   1.6999993      -133.53346    
   1.7099993      0.21409836    
   1.7099993      -8.1679096    
   1.7099993      -133.89883    
   1.7199993      0.23839833    
   1.7199993      -8.0294447    
   1.7199993      -134.18831    
   1.7299993      0.26289830    
   1.7299993      -7.8863430    
   1.7299993      -134.39798    
   1.7399993      0.28759825    
   1.7399993      -7.7385349    
   1.7399993      -134.52386    
   1.7499993      0.31249821    
   1.7499993      -7.5859485    
   1.7499993      -134.56177    
   1.7599993      0.33759817    
   1.7599993      -7.4285121    
   1.7599993      -134.50745    
   1.7699993      0.36289814    
   1.7699993      -7.2661533    
   1.7699993      -134.35651    
   1.7799993      0.38839808    
   1.7799993      -7.0987997    
   1.7799993      -134.10443    
   1.7899992      0.41409805    
   1.7899992      -6.9263778    
   1.7899992      -133.74652    
   1.7999992      0.43999803    
   1.7999992      -6.7488136    
   1.7999992      -133.27805    
   1.8099992      0.46609798    
   1.8099992      -6.5660338    
   1.8099992      -132.69403    
   1.8199992      0.49239793    
   1.8199992      -6.3779635    
   1.8199992      -131.98946    
   1.8299992      0.51889789    
   1.8299992      -6.1845269    
   1.8299992      -131.15907    
   1.8399992      0.54559785    
   1.8399992      -5.9856496    
   1.8399992      -130.19754    
   1.8499992      0.57249779    
   1.8499992      -5.7812543    
   1.8499992      -129.09935    
   1.8599992      0.59959775    
   1.8599992      -5.5712652    
   1.8599992      -127.85886    
   1.8699992      0.62689775    
   1.8699992      -5.3556046    
   1.8699992      -126.47026    
   1.8799992      0.65439767    
   1.8799992      -5.1341963    
   1.8799992      -124.92759    
   1.8899992      0.68209761    
   1.8899992      -4.9069610    
   1.8899992      -123.22472    
   1.8999991      0.70999759    
   1.8999991      -4.6738205    
   1.8999991      -121.35535    
   1.9099991      0.73809755    
   1.9099991      -4.4346957    
   1.9099991      -119.31307    
   1.9199991      0.76639754    
   1.9199991      -4.1895075    
   1.9199991      -117.09122    
   1.9299991      0.79489744    
   1.9299991      -3.9381766    
   1.9299991      -114.68305    
   1.9399991      0.82359743    
   1.9399991      -3.6806214    
   1.9399991      -112.08156    
   1.9499991      0.85249740    
   1.9499991      -3.4167616    
   1.9499991      -109.27963    
   1.9599991      0.88159734    
   1.9599991      -3.1465158    
   1.9599991      -106.26992    
   1.9699991      0.91089725    
   1.9699991      -2.8698025    
   1.9699991      -103.04493    
   1.9799991      0.94039726    
   1.9799991      -2.5865383    
   1.9799991      -99.596970    
   1.9899991      0.97009718    
   1.9899991      -2.2966418    
   1.9899991      -95.918152    
   1.9999990      0.99999714    
   1.9999990      -2.0000286    
   1.9999990      -92.000381    
   2.0099990       1.0300971    
   2.0099990      -1.6966155    
   2.0099990      -87.835403    
   2.0199990       1.0603970    
   2.0199990      -1.3863186    
   2.0199990      -83.414749    
   2.0299990       1.0908970    
   2.0299990      -1.0690522    
   2.0299990      -78.729706    
   2.0399990       1.1215969    
   2.0399990     -0.74473161    
   2.0399990      -73.771393    
   2.0499990       1.1524969    
   2.0499990     -0.41327053    
   2.0499990      -68.530716    
   2.0599990       1.1835968    
   2.0599990     -7.45845437E-02
   2.0599990      -62.998379    
   2.0699990       1.2148968    
   2.0699990      0.27141431    
   2.0699990      -57.164829    
   2.0799990       1.2463968    
   2.0799990      0.62481344    
   2.0799990      -51.020306    
   2.0899990       1.2780967    
   2.0899990      0.98569936    
   2.0899990      -44.554848    
   2.0999990       1.3099966    
   2.0999990       1.3541604    
   2.0999990      -37.758247    
   2.1099989       1.3420966    
   2.1099989       1.7302866    
   2.1099989      -30.620033    
   2.1199989       1.3743966    
   2.1199989       2.1141658    
   2.1199989      -23.129545    
   2.1299989       1.4068965    
   2.1299989       2.5058866    
   2.1299989      -15.275887    
   2.1399989       1.4395964    
   2.1399989       2.9055402    
   2.1399989      -7.0478663    
   2.1499989       1.4724964    
   2.1499989       3.3132176    
   2.1499989       1.5659319    
   2.1599989       1.5055963    
   2.1599989       3.7290082    
   2.1599989       10.577140    
   2.1699989       1.5388963    
   2.1699989       4.1530056    
   2.1699989       19.997725    
   2.1799989       1.5723963    
   2.1799989       4.5852990    
   2.1799989       29.839811    
   2.1899989       1.6060961    
   2.1899989       5.0259819    
   2.1899989       40.115864    
   2.1999989       1.6399961    
   2.1999989       5.4751482    
   2.1999989       50.838604    
   2.2099988       1.6740961    
   2.2099988       5.9328904    
   2.2099988       62.021027    
   2.2199988       1.7083960    
   2.2199988       6.3993025    
   2.2199988       73.676361    
   2.2299988       1.7428960    
   2.2299988       6.8744793    
   2.2299988       85.818192    
   2.2399988       1.7775959    
   2.2399988       7.3585138    
   2.2399988       98.460289    
   2.2499988       1.8124958    
   2.2499988       7.8515038    
   2.2499988       111.61683    
   2.2599988       1.8475958    
   2.2599988       8.3535433    
   2.2599988       125.30217    
   2.2699988       1.8828957

Notice that I entered h=0.01.
My error message is

*** glibc detected *** ./a.out: free(): invalid next size (fast): 0x091e5018 ***
======= Backtrace: =========
/lib/i386-linux-gnu/libc.so.6(+0x6b961)[0x236961]
/lib/i386-linux-gnu/libc.so.6(+0x6d28b)[0x23828b]
/lib/i386-linux-gnu/libc.so.6(cfree+0x6d)[0x23b41d]
./a.out[0x804901d]
./a.out[0x8049093]
/lib/i386-linux-gnu/libc.so.6(__libc_start_main+0xe7)[0x1e1e37]
./a.out[0x8048761]
======= Memory map: ========
00110000-001c8000 r-xp 00000000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
001c8000-001c9000 r--p 000b7000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
001c9000-001ca000 rw-p 000b8000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
001ca000-001cb000 rw-p 00000000 00:00 0 
001cb000-00325000 r-xp 00000000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00325000-00326000 ---p 0015a000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00326000-00328000 r--p 0015a000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00328000-00329000 rw-p 0015c000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00329000-0032c000 rw-p 00000000 00:00 0 
0061e000-0061f000 r-xp 00000000 00:00 0          [vdso]
00ae0000-00b04000 r-xp 00000000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00b04000-00b05000 r--p 00023000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00b05000-00b06000 rw-p 00024000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00c5f000-00c79000 r-xp 00000000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00c79000-00c7a000 r--p 00019000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00c7a000-00c7b000 rw-p 0001a000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00ea9000-00ec5000 r-xp 00000000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00ec5000-00ec6000 r--p 0001b000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00ec6000-00ec7000 rw-p 0001c000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
08048000-0804a000 r-xp 00000000 08:06 529439     /home/fluidistic/numerical2011/guide5/a.out
0804a000-0804b000 r--p 00001000 08:06 529439     /home/fluidistic/numerical2011/guide5/a.out
0804b000-0804c000 rw-p 00002000 08:06 529439     /home/fluidistic/numerical2011/guide5/a.out
0804c000-08064000 rw-p 00000000 00:00 0 
091c9000-091ea000 rw-p 00000000 00:00 0          [heap]
b7700000-b7721000 rw-p 00000000 00:00 0 
b7721000-b7800000 ---p 00000000 00:00 0 
b78ce000-b78d0000 rw-p 00000000 00:00 0 
b78ed000-b78ef000 rw-p 00000000 00:00 0 
bfc6b000-bfc8c000 rw-p 00000000 00:00 0          [stack]
Aborted

 

[Mark44]

Thanks for your help.
I tried to define new variables (scalars) in function of x(i), y(i) in the subroutine as to get rid of the arrays, but without success.

That's the wrong direction to go. Instead of changing x and y to scalars in your main program, change the declaration of d in your Newton subroutine so that it is an array.

These errors are telling you that there is something wrong with how d is declared in the subroutine.

Newton.f90:34.2:

  d(i)=y(i)
  1
Error: Unclassifiable statement at (1)
Newton.f90:38.12:

            d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
            1
Error: Unclassifiable statement at (1)
Newton.f90:24.27:

subroutine Newton(x_m,h,x,d,y)
                           1
Error: PROCEDURE attribute conflicts with INTENT attribute in 'd' at (1)
Newton.f90:41.11:

I mistakenly said that x was being used as a scalar in your subroutine, but that is not the case. However, I stand by what I said about your d variable.

Think about how you're declaring d in the main program, and how it shows up when you call your subroutine, and then how you have the formal variable of the same name in the subroutine declared: (I've chopped out many lines of code.)

real, dimension(:), allocatable :: x,y,d,p
real :: h, x_m
.
.
allocate(x(0:n),y(0:n), d(1:n), p(1:n)) // x, y, d, and p are arrays
.
.

call Newton(x_m,h,x,d,y) // Actual args are scalar, scalar, array, array, array
.
.
// In line below, the formal parameters are scalar, scalar, array, [color="red"]scalar[/color], array. 
// That is what you need to fix -- make d an array.
subroutine Newton(x_m,h,x,d,y) 
real, intent(inout)::x_m,h,d
real:: prod, p
real, dimension(:), intent(inout) :: x,y

I'm really confused about how to solve the problem. Basically I must transform the arrays x and d into scalars, in the subroutine?

Edit: I tried to modify only a small part of the declaration of variables in the subroutine.
Now it looks like this:

subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d

The program compiles! But... I do get an error.
Here it is:

Enter the degree of the polynomial n and a step-size h
3, 0.001
 Enter the name of the file to read the data from
kincaid.dat
At line 10 of file hope2.f90 (unit = 5, file = 'stdin')
Fortran runtime error: Insufficient data descriptors in format after reversion
fluidistic@fluidistic-desktop:~/numerical2011/guide5$ ncaid.dat
/usr/bin/python: can't find '__main__' module in '/usr/share/command-not-found'

My line 10 is

read(*,'(x)')datos

I'm guessing the '(x)' part is a problem. I must say I randomly wrote '(x)'. My professor gave us a -working- program with such an example.
It starts like this:

character(50)::matrix, vector
  real, dimension(:,:), allocatable::A
  real, dimension(:), allocatable::b,x,order
  integer::n,i,j

  write(*,*)'Enter n'
  read(*,*)n

  write(*,*)'Enter the name of the file of the matrix A'
  read(*,'(A)')matrix
  open(10,file=matrix,status='old')

  write(*,*)'Do the same for the b vector'
  read(*,'(A)')vector
  open(15,file=vector,status='old')

  allocate(A(n,n),b(n),x(n),order(n))

But I really do not understand how the read(*,'whatisthis?') works.

Edit 2: I searched on the internet and I've learned that '(A)' has nothing to do with the A matrix. Instead, it's a format or something like that. So I replaced my '(x)' by '(A)'. I also erased some "".
Now my program looks like this:

program div
implicit none
character(100000) :: datos
real, dimension(:), allocatable :: x,y,d,p
integer::n ,i
real :: h, x_m
write(*,*)"Enter the degree of the polynomial n and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(unit=1,file=datos, status="old")
allocate(x(0:n),y(0:n), d(1:n), p(1:n))
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)

call Newton(x_m,h,x,d,y)

contains
subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=120, file="polynomial.dat")
do while (x_m <= x(n))
  do i=0,n
     d(i)=y(i)
  end do
  do j=1, n
     if (j>= i) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
  prod=x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*x_m-x(j)
    end do
    p=d(i)*prod
    write(120,*)p
  end do    

  x_m=x_m+h
end do
close(120)
end subroutine
end program

But I'm having an error of the form *** glibc detected *** ./a.out: free(): invalid next size (fast): 0x08f74018 ***.

I'm guessing it has to do with the '(A)' format?
I'm trying to read from a file with 2 columns of n rows. What should I use instead of '(A)'?

Edit3: I could compile+execute without error, though I had to modificate badly my program. I get only 1 number as output instead of more than 100.
I'll work on that...

 

[fluidistic]

That's the wrong direction to go. Instead of changing x and y to scalars in your main program, change the declaration of d in your Newton subroutine so that it is an array.

These errors are telling you that there is something wrong with how d is declared in the subroutine.

Thank you once again for your help.
I know my last 2 posts are a total mess. I think it's safe if you only read my last post. I fixed the error of the arrays/scalars already.
Meanwhile I just see that I wrongly placed my "write" statement for x_m and p. I'll try to fix this now.
Feel free to have a look over my last post.

Edit: I moved the write statement outside the loop it was in. Now I get decent values for x_m but still wrong values for p. And I still get the memory error message. I wonder what's wrong. It likely has to see with the files I open/close.

 

[Mark44]

Here are some things I see:
1) character(100000) :: datos
This is saying that the name of this file will be 100,000 characters long.
2)The 2nd and 3rd line below are in the wrong order.

write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(unit=1,file=datos, status="old")

First, you open the file. Next, you read from it.
The 1st parameter in the read statement is the unit number (1), and the second parameter can be the line number of a format statement.
Also, datos is supposed to be the name of your file, not a variable that you will store a number in.

Do a web search for "fortran read format" -- that should probably give you a bunch of links to info on how to do formatted input from a file.

 

[Mark44]

One other thing. Since you seem to be getting bogus numbers (at least in post #4), try your program with a much larger step size, say 0.5 rather than .001. Then you can verify that your program is producing reasonable numbers.

 

[fluidistic]

Here are some things I see:
1) character(100000) :: datos
This is saying that the name of this file will be 100,000 characters long.

I changed it to 5 since datos is 5 leters long but the program would crash because it can't read "kincaid.dat". Setting character (50) or character (100000) is the same for me. My program can actually read "kincaid.dat" with character 1000000 :) I find it safer to put a large value, just in case.

2)The 2nd and 3rd line below are in the wrong order.

write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(unit=1,file=datos, status="old")

First, you open the file. Next, you read from it.

Actually no. Let me explain.
What happens with

write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(1,file=datos)

is: I input "kincaid.dat". My program call it datos, I mean it assignates kincaid.dat as datos. So the line read(*,'(A)')datos actually reads the file kincaid.dat, not datos since there's no such a file.
Then the line open(1,file=datos) does not create the file "datos" nor does it open the file "datos" since there's no such a file. Instead, it opens the file "kincaid.dat".
I do know this works because my x(0) is indeed the one of kincaid.dat.
I even tried to invert these lines and I got an error like "datos" doesn't exist.

The 1st parameter in the read statement is the unit number (1), and the second parameter can be the line number of a format statement.
Also, datos is supposed to be the name of your file, not a variable that you will store a number in.

Do a web search for "fortran read format" -- that should probably give you a bunch of links to info on how to do formatted input from a file.

I did read about "read format". I've learned that when dealing with a character variable, I must use '(A)'. I could specify 2 other options about it, but I don't think they are relevant. (like reading the first 3 digits of each line in "kincaid.dat", etc.)

One other thing. Since you seem to be getting bogus numbers (at least in post #4), try your program with a much larger step size, say 0.5 rather than .001. Then you can verify that your program is producing reasonable numbers.

Well I moved the write x_m, p out of its do loop.
Now I get as "polynomial.dat"

   1.0000000      -56.000000    
   1.0100000      -56.888004    
   1.0200000      -57.792053    
   1.0300000      -58.712166    
   1.0400000      -59.648361    
   1.0500000      -60.600628    
   1.0599999      -61.568951    
   1.0699999      -62.553299    
   1.0799999      -63.553600    
   1.0899999      -64.569801    
   1.0999999      -65.601807    
   1.1099999      -66.649490    
   1.1199999      -67.712730    
   1.1299999      -68.791367    
   1.1399999      -69.885223    
   1.1499999      -70.994102    
   1.1599998      -72.117767    
   1.1699998      -73.255966    
   1.1799998      -74.408432    
   1.1899998      -75.574837    
   1.1999998      -76.754868    
   1.2099998      -77.948151    
   1.2199998      -79.154289    
   1.2299998      -80.372864    
   1.2399998      -81.603409    
   1.2499998      -82.845428    
   1.2599998      -84.098404    
   1.2699997      -85.361755    
   1.2799997      -86.634895    
   1.2899997      -87.917175    
   1.2999997      -89.207909    
   1.3099997      -90.506401    
   1.3199997      -91.811852    
   1.3299997      -93.123482    
   1.3399997      -94.440430    
   1.3499997      -95.761787    
   1.3599997      -97.086632    
   1.3699996      -98.413948    
   1.3799996      -99.742691    
   1.3899996      -101.07178    
   1.3999996      -102.40005    
   1.4099996      -103.72630    
   1.4199996      -105.04927    
   1.4299996      -106.36764    
   1.4399996      -107.68003    
   1.4499996      -108.98500    
   1.4599996      -110.28105    
   1.4699996      -111.56660    
   1.4799995      -112.84005    
   1.4899995      -114.09967    
   1.4999995      -115.34370    
   1.5099995      -116.57029    
   1.5199995      -117.77753    
   1.5299995      -118.96346    
   1.5399995      -120.12599    
   1.5499995      -121.26297    
   1.5599995      -122.37219    
   1.5699995      -123.45133    
   1.5799994      -124.49804    
   1.5899994      -125.50980    
   1.5999994      -126.48407    
   1.6099994      -127.41818    
   1.6199994      -128.30942    
   1.6299994      -129.15489    
   1.6399994      -129.95172    
   1.6499994      -130.69684    
   1.6599994      -131.38710    
   1.6699994      -132.01930    
   1.6799994      -132.59009    
   1.6899993      -133.09599    
   1.6999993      -133.53346    
   1.7099993      -133.89883    
   1.7199993      -134.18831    
   1.7299993      -134.39798    
   1.7399993      -134.52386    
   1.7499993      -134.56177    
   1.7599993      -134.50745    
   1.7699993      -134.35651    
   1.7799993      -134.10443    
   1.7899992      -133.74652    
   1.7999992      -133.27805    
   1.8099992      -132.69403    
   1.8199992      -131.98946    
   1.8299992      -131.15907    
   1.8399992      -130.19754    
   1.8499992      -129.09935    
   1.8599992      -127.85886    
   1.8699992      -126.47026    
   1.8799992      -124.92759    
   1.8899992      -123.22472    
   1.8999991      -121.35535    
   1.9099991      -119.31307    
   1.9199991      -117.09122    
   1.9299991      -114.68305    
   1.9399991      -112.08156    
   1.9499991      -109.27963    
   1.9599991      -106.26992    
   1.9699991      -103.04493    
   1.9799991      -99.596970    
   1.9899991      -95.918152    
   1.9999990      -92.000381    
   2.0099990      -87.835403    
   2.0199990      -83.414749    
   2.0299990      -78.729706    
   2.0399990      -73.771393    
   2.0499990      -68.530716    
   2.0599990      -62.998379    
   2.0699990      -57.164829    
   2.0799990      -51.020306    
   2.0899990      -44.554848    
   2.0999990      -37.758247    
   2.1099989      -30.620033    
   2.1199989      -23.129545    
   2.1299989      -15.275887    
   2.1399989      -7.0478663    
   2.1499989       1.5659319    
   2.1599989       10.577140    
   2.1699989       19.997725    
   2.1799989       29.839811    
   2.1899989       40.115864    
   2.1999989       50.838604    
   2.2099988       62.021027    
   2.2199988       73.676361    
   2.2299988       85.818192    
   2.2399988       98.460289    
   2.2499988       111.61683    
   2.2599988       125.30217    
   2.2699988       139.53101    
   2.2799988       154.31844    
   2.2899988       169.67975    
   2.2999988       185.63046    
   2.3099988       202.18666    
   2.3199987       219.36456    
   2.3299987       237.18074    
   2.3399987       255.65216    
   2.3499987       274.79599    
   2.3599987       294.62997    
   2.3699987       315.17181    
   2.3799987       336.44000    
   2.3899987       358.45313    
   2.3999987       381.23022    
   2.4099987       404.79062    
   2.4199986       429.15393    
   2.4299986       454.34052    
   2.4399986       480.37061    
   2.4499986       507.26529    
   2.4599986       535.04565    
   2.4699986       563.73352    
   2.4799986       593.35083    
   2.4899986       623.92010    
   2.4999986       655.46417    
   2.5099986       688.00635    
   2.5199986       721.57025    
   2.5299985       756.18036    
   2.5399985       791.86066    
   2.5499985       828.63684    
   2.5599985       866.53387    
   2.5699985       905.57782    
   2.5799985       945.79498    
   2.5899985       987.21228    
   2.5999985       1029.8569    
   2.6099985       1073.7570    
   2.6199985       1118.9399    
   2.6299984       1165.4354    
   2.6399984       1213.2722    
   2.6499984       1262.4801    
   2.6599984       1313.0891    
   2.6699984       1365.1305    
   2.6799984       1418.6349    
   2.6899984       1473.6344    
   2.6999984       1530.1613    
   2.7099984       1588.2485    
   2.7199984       1647.9294    
   2.7299984       1709.2377    
   2.7399983       1772.2081    
   2.7499983       1836.8756    
   2.7599983       1903.2758    
   2.7699983       1971.4451    
   2.7799983       2041.4194    
   2.7899983       2113.2375    
   2.7999983       2186.9365    
   2.8099983       2262.5547    
   2.8199983       2340.1321    
   2.8299983       2419.7075    
   2.8399982       2501.3228    
   2.8499982       2585.0173    
   2.8599982       2670.8337    
   2.8699982       2758.8142    
   2.8799982       2849.0020    
   2.8899982       2941.4402    
   2.8999982       3036.1736    
   2.9099982       3133.2466    
   2.9199982       3232.7058    
   2.9299982       3334.5972    
   2.9399981       3438.9670    
   2.9499981       3545.8643    
   2.9599981       3655.3369    
   2.9699981       3767.4343    
   2.9799981       3882.2063    
   2.9899981       3999.7031    
   2.9999981       4119.9771    
   3.0099981       4243.0791    
   3.0199981       4369.0640    
   3.0299981       4497.9834    
   3.0399981       4629.8931    
   3.0499980       4764.8491    
   3.0599980       4902.9048    
   3.0699980       5044.1196    
   3.0799980       5188.5498    
   3.0899980       5336.2549    
   3.0999980       5487.2939    
   3.1099980       5641.7261    
   3.1199980       5799.6143    
   3.1299980       5961.0190    
   3.1399980       6126.0034    
   3.1499979       6294.6304    
   3.1599979       6466.9658    
   3.1699979       6643.0742    
   3.1799979       6823.0210    
   3.1899979       7006.8760    
   3.1999979       7194.7046    
   3.2099979       7386.5771    
   3.2199979       7582.5645    
   3.2299979       7782.7349    
   3.2399979       7987.1621    
   3.2499979       8195.9199    
   3.2599978       8409.0801    
   3.2699978       8626.7178    
   3.2799978       8848.9111    
   3.2899978       9075.7344    
   3.2999978       9307.2656    
   3.3099978       9543.5840    
   3.3199978       9784.7705    
   3.3299978       10030.905    
   3.3399978       10282.070    
   3.3499978       10538.348    
   3.3599977       10799.823    
   3.3699977       11066.580    
   3.3799977       11338.707    
   3.3899977       11616.290    
   3.3999977       11899.417    
   3.4099977       12188.180    
   3.4199977       12482.667    
   3.4299977       12782.972    
   3.4399977       13089.187    
   3.4499977       13401.403    
   3.4599977       13719.723    
   3.4699976       14044.239    
   3.4799976       14375.049    
   3.4899976       14712.251    
   3.4999976       15055.947    
   3.5099976       15406.240    
   3.5199976       15763.230    
   3.5299976       16127.022    
   3.5399976       16497.723    
   3.5499976       16875.436    
   3.5599976       17260.270    
   3.5699975       17652.336    
   3.5799975       18051.742    
   3.5899975       18458.600    
   3.5999975       18873.023    
   3.6099975       19295.135    
   3.6199975       19725.037    
   3.6299975       20162.854    
   3.6399975       20608.701    
   3.6499975       21062.707    
   3.6599975       21524.980    
   3.6699975       21995.650    
   3.6799974       22474.846    
   3.6899974       22962.684    
   3.6999974       23459.299    
   3.7099974       23964.813    
   3.7199974       24479.363    
   3.7299974       25003.076    
   3.7399974       25536.086    
   3.7499974       26078.527    
   3.7599974       26630.537    
   3.7699974       27192.256    
   3.7799973       27763.818    
   3.7899973       28345.367    
   3.7999973       28937.043    
   3.8099973       29538.996    
   3.8199973       30151.363    
   3.8299973       30774.299    
   3.8399973       31407.947    
   3.8499973       32052.463    
   3.8599973       32707.998    
   3.8699973       33374.703    
   3.8799973       34052.738    
   3.8899972       34742.258    
   3.8999972       35443.422    
   3.9099972       36156.391    
   3.9199972       36881.324    
   3.9299972       37618.395    
   3.9399972       38367.758    
   3.9499972       39129.590    
   3.9599972       39904.063    
   3.9699972       40691.340    
   3.9799972       41491.598    
   3.9899971       42305.008    
   3.9999971       43131.762    
   4.0099974       43972.043    
   4.0199976       44826.020    
   4.0299978       45693.879    
   4.0399981       46575.789    
   4.0499983       47471.953    
   4.0599985       48382.559    
   4.0699987       49307.789    
   4.0799990       50247.844    
   4.0899992       51202.910    
   4.0999994       52173.191    
   4.1099997       53158.891    
   4.1199999       54160.188    
   4.1300001       55177.309    
   4.1400003       56210.445    
   4.1500006       57259.816    
   4.1600008       58325.625    
   4.1700010       59408.074    
   4.1800013       60507.387    
   4.1900015       61623.773    
   4.2000017       62757.465    
   4.2100019       63908.652    
   4.2200022       65077.594    
   4.2300024       66264.492    
   4.2400026       67469.578    
   4.2500029       68693.078    
   4.2600031       69935.227    
   4.2700033       71196.266    
   4.2800035       72476.414    
   4.2900038       73775.922    
   4.3000040       75095.023    
   4.3100042       76433.961    
   4.3200045       77792.992    
   4.3300047       79172.352    
   4.3400049       80572.289    
   4.3500051       81993.070    
   4.3600054       83434.938    
   4.3700056       84898.156    
   4.3800058       86382.969    
   4.3900061       87889.664    
   4.4000063       89418.500    
   4.4100065       90969.734    
   4.4200068       92543.641    
   4.4300070       94140.492    
   4.4400072       95760.555    
   4.4500074       97404.133    
   4.4600077       99071.484    
   4.4700079       100762.89    
   4.4800081       102478.64    
   4.4900084       104219.05    
   4.5000086       105984.38    
   4.5100088       107774.93    
   4.5200090       109590.98    
   4.5300093       111432.87    
   4.5400095       113300.88    
   4.5500097       115195.30    
   4.5600100       117116.47    
   4.5700102       119064.68    
   4.5800104       121040.24    
   4.5900106       123043.49    
   4.6000109       125074.73    
   4.6100111       127134.27    
   4.6200113       129222.47    
   4.6300116       131339.63    
   4.6400118       133486.11    
   4.6500120       135662.20    
   4.6600122       137868.28    
   4.6700125       140104.67    
   4.6800127       142371.72    
   4.6900129       144669.77    
   4.7000132       146999.14    
   4.7100134       149360.27    
   4.7200136       151753.45    
   4.7300138       154179.03    
   4.7400141       156637.39    
   4.7500143       159128.91    
   4.7600145       161653.94    
   4.7700148       164212.88    
   4.7800150       166806.05    
   4.7900152       169433.86    
   4.8000154       172096.73    
   4.8100157       174794.98    
   4.8200159       177529.03    
   4.8300161       180299.30    
   4.8400164       183106.11    
   4.8500166       185949.95    
   4.8600168       188831.13    
   4.8700171       191750.13    
   4.8800173       194707.30    
   4.8900175       197703.11    
   4.9000177       200737.95    
   4.9100180       203812.23    
   4.9200182       206926.36    
   4.9300184       210080.83    
   4.9400187       213276.00    
   4.9500189       216512.31    
   4.9600191       219790.28    
   4.9700193       223110.27    
   4.9800196       226472.75    
   4.9900198       229878.17    
   5.0000200       233326.97    
   5.0100203       236819.59    
   5.0200205       240356.55    
   5.0300207       243938.27    
   5.0400209       247565.22    
   5.0500212       251237.89    
   5.0600214       254956.77    
   5.0700216       258722.30    
   5.0800219       262535.00    
   5.0900221       266395.28    
   5.1000223       270303.78    
   5.1100225       274260.81    
   5.1200228       278267.03    
   5.1300230       282322.84    
   5.1400232       286428.78    
   5.1500235       290585.47    
   5.1600237       294793.22    
   5.1700239       299052.69    
   5.1800241       303364.41    
   5.1900244       307728.78    
   5.2000246       312146.50    
   5.2100248       316618.03    
   5.2200251       321143.88    
   5.2300253       325724.63    
   5.2400255       330360.88    
   5.2500257       335053.06    
   5.2600260       339801.88    
   5.2700262       344607.78    
   5.2800264       349471.38    
   5.2900267       354393.28    
   5.3000269       359374.00    
   5.3100271       364414.13    
   5.3200274       369514.28    
   5.3300276       374675.06    
   5.3400278       379897.03    
   5.3500280       385180.81    
   5.3600283       390527.00    
   5.3700285       395936.19    
   5.3800287       401409.03    
   5.3900290       406946.09    
   5.4000292       412548.06    
   5.4100294       418215.56    
   5.4200296       423949.06    
   5.4300299       429749.47    
   5.4400301       435617.22    
   5.4500303       441553.03    
   5.4600306       447557.66    
   5.4700308       453631.56    
   5.4800310       459775.53    
   5.4900312       465990.19    
   5.5000315       472276.28    
   5.5100317       478634.47    
   5.5200319       485065.31    
   5.5300322       491569.66    
   5.5400324       498148.06    
   5.5500326       504801.25    
   5.5600328       511530.13    
   5.5700331       518335.13    
   5.5800333       525217.19    
   5.5900335       532176.81    
   5.6000338       539214.75    
   5.6100340       546332.00    
   5.6200342       553529.06    
   5.6300344       560806.75    
   5.6400347       568165.75    
   5.6500349       575606.88    
   5.6600351       583130.94    
   5.6700354       590738.50    
   5.6800356       598430.56    
   5.6900358       606207.81    
   5.7000360       614071.00    
   5.7100363       622020.94    
   5.7200365       630058.31    
   5.7300367       638184.13    
   5.7400370       646399.06    
   5.7500372       654703.94    
   5.7600374       663099.56    
   5.7700377       671586.81    
   5.7800379       680166.38    
   5.7900381       688839.31    
   5.8000383       697606.38    
   5.8100386       706468.19    
   5.8200388       715425.94    
   5.8300390       724480.25    
   5.8400393       733632.00    
   5.8500395       742882.25    
   5.8600397       752231.81    
   5.8700399       761681.38    
   5.8800402       771231.88    
   5.8900404       780884.50    
   5.9000406       790639.75    
   5.9100409       800498.94    
   5.9200411       810462.56    
   5.9300413       820531.75    
   5.9400415       830707.50    
   5.9500418       840990.69    
   5.9600420       851382.25    
   5.9700422       861883.00    
   5.9800425       872494.00    
   5.9900427       883216.31

So basically my x_m's are right and my p's are wrong. I took the algorithm from Kincaid's book... I really don't know what the heck is wrong with it.
And what about this memory error I get? :/
Seems like my program compiles perfectly, runs entirely but create a memory error and my output values p are all wrong.Edit: I modified my subroutine. I still get non sense result for p+memory error but I'm getting closer to a correct program.
I have

subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
 do i=0,n
    d(i)=y(i)
! end do
  do j=1, n
     if (i>= j) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
end do
  prod=1.!x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*(x_m-x(j))
    end do
    p=d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine
end program

 

[Mark44]

I changed it to 5 since datos is 5 leters long but the program would crash because it can't read "kincaid.dat".

Of course not. The string "kincaid.dat" has 11 characters, so you can't store it in an array that holds only 5. You are completely misunderstanding the difference between a variable and what is stored in it. You are not going to store the name of the variable in the variable - you are going to store the name of the file. Since the file name you are using has 11 characters, counting the period, your variable needs to be at least that long. Declaring a character array variable that can hold 50 characters is reasonable - declaring one to hold 100,000 characters is not reasonable.

Setting character (50) or character (100000) is the same for me. My program can actually read "kincaid.dat" with character 1000000 :) I find it safer to put a large value, just in case.

Actually no. Let me explain.
What happens with

write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(1,file=datos)

is: I input "kincaid.dat". My program call it datos, I mean it assignates kincaid.dat as datos. So the line read(*,'(A)')datos actually reads the file kincaid.dat, not datos since there's no such a file.
Then the line open(1,file=datos) does not create the file "datos" nor does it open the file "datos" since there's no such a file. Instead, it opens the file "kincaid.dat".
I do know this works because my x(0) is indeed the one of kincaid.dat.
I even tried to invert these lines and I got an error like "datos" doesn't exist.

I misunderstood what you were trying to do before, so I take back what I said about the lines being reversed.

Your explanation above is off, though.
The prompt (the write statement) tells the user to enter the name of the file.
The read statement reads some characters from the keyboard and stores them in the datos character array. It DOES NOT read the file.
The open statement associates unit # 1 with the file whose name is what you typed in for the read statement.

I did read about "read format". I've learned that when dealing with a character variable, I must use '(A)'. I could specify 2 other options about it, but I don't think they are relevant. (like reading the first 3 digits of each line in "kincaid.dat", etc.)

Well I moved the write x_m, p out of its do loop.
Now I get as "polynomial.dat"

   1.0000000      -56.000000    
   1.0100000      -56.888004    
   1.0200000      -57.792053    
   1.0300000      -58.712166       
   <cut>   
   5.9800425       872494.00    
   5.9900427       883216.31

So basically my x_m's are right and my p's are wrong.

I really didn't need to see 500 lines of output that may or may not be garbage.

I took the algorithm from Kincaid's book... I really don't know what the heck is wrong with it.
And what about this memory error I get? :/
Seems like my program compiles perfectly, runs entirely but create a memory error and my output values p are all wrong.

I don't know Kincaid's book and you've shown only a few lines of the algorithm you're using. If you're going to implement an algorithm in some language (like Fortran) you need to understand what is being said in the algorithm, so let's get back to that. Please include the full algorithm so we can see if your code is an accurate rendition of it.

What memory error are you talking about? You have mentioned it, but each of your posts is fairly long, and I don't think you gave much of a description about it. In any case, you need to make sure that you have the algorithm right.

Edit: I modified my subroutine. I still get non sense result for p+memory error but I'm getting closer to a correct program.
I have

subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
 do i=0,n
    d(i)=y(i)
! end do
  do j=1, n
     if (i>= j) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
end do
  prod=1.!x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*(x_m-x(j))
    end do
    p=d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine
end program

 

[fluidistic]

Thanks again for your help!

Of course not. The string "kincaid.dat" has 11 characters, so you can't store it in an array that holds only 5. You are completely misunderstanding the difference between a variable and what is stored in it. You are not going to store the name of the variable in the variable - you are going to store the name of the file. Since the file name you are using has 11 characters, counting the period, your variable needs to be at least that long. Declaring a character array variable that can hold 50 characters is reasonable - declaring one to hold 100,000 characters is not reasonable.

Ok I understand well now. I've set it to 50 in case I want to use this program with another file than kincaid.dat.

I misunderstood what you were trying to do before, so I take back what I said about the lines being reversed.

Your explanation above is off, though.
The prompt (the write statement) tells the user to enter the name of the file.
The read statement reads some characters from the keyboard and stores them in the datos character array. It DOES NOT read the file.
The open statement associates unit # 1 with the file whose name is what you typed in for the read statement.

Nice, I didn't realize.

I really didn't need to see 500 lines of output that may or may not be garbage.

You're right. My next will be 6 lines.

I don't know Kincaid's book and you've shown only a few lines of the algorithm you're using. If you're going to implement an algorithm in some language (like Fortran) you need to understand what is being said in the algorithm, so let's get back to that. Please include the full algorithm so we can see if your code is an accurate rendition of it.

Roger that.

What memory error are you talking about? You have mentioned it, but each of your posts is fairly long, and I don't think you gave much of a description about it. In any case, you need to make sure that you have the algorithm right.

Ok the algorithm given isn't in Fortran.
What Kincaid does is first calculate the coefficient that he calls "di" such that the interpolating polynomial is worth $p(x)= \sum _{i=0}^n di \prod _{j=0}^{i-1} (x-x_j)$.
He says that the most efficient way to get those di is to use the following algorithm:

for i=0, 1,...,n do
di <- f(xi)
end
for j=1,2,...,n do
for i=j, j+1, ..., n do
di <-(di-d_{i-1} )/(xi-x_{i-j})
end 
end.

This is the entire algorithm. For a fixed di, I do not understand exactly why the next coefficient will be "(di-d_{i-1} )/(xi-x_{i-j})" but I think it has to see with divided differences.
I understand why and that the first di=f(x0).
In my first posts I set this algorithm almost without thinking, hence my wrong values for p. But now if you look at my current code, I do not understand why I'm getting wrong values.
By the way I agree with you that my posts aren't clear. I'll post the latest updates now and here.
My code:

program div
implicit none
character(50) :: datos
real, dimension(:), allocatable :: x,y,d
integer::n ,i
real :: h, x_m
write(*,*)"Enter the degree of the polynomial n and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file to read the data from"
read(*,'(A)')datos
open(1,file=datos)

allocate(x(0:n),y(0:n), d(1:n))
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)
x_m=x(0)
call Newton(x_m,h,x,d,y)

contains
subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
 do i=0,n
    d(i)=y(i)
! end do
  do j=1, n
     if (i>= j) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
end do
  prod=1.!x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*(x_m-x(j))
    end do
    p=d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine
end program

My kincaid.dat file:

1 -3
3 1
5 2
6 4

I'll now show you my input+error message I get. (I copy and paste from the terminal in which I compile + execute the program).

fluidistic@fluidistic-desktop:~/numerical2011/guide5$ ./a.out
 Enter the degree of the polynomial n and a step-size h
3, 1
 Enter the name of the file to read the data from
kincaid.dat
*** glibc detected *** ./a.out: free(): invalid next size (fast): 0x093799b0 ***
======= Backtrace: =========
/lib/i386-linux-gnu/libc.so.6(+0x6b961)[0x4a5961]
/lib/i386-linux-gnu/libc.so.6(+0x6d28b)[0x4a728b]
/lib/i386-linux-gnu/libc.so.6(cfree+0x6d)[0x4aa41d]                                                                                                                     
./a.out[0x8049057]
./a.out[0x80490d0]                                                                                                                                                      
/lib/i386-linux-gnu/libc.so.6(__libc_start_main+0xe7)[0x450e37]
./a.out[0x8048761]
======= Memory map: ========                                                                                                                                            
0019a000-00252000 r-xp 00000000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
00252000-00253000 r--p 000b7000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
00253000-00254000 rw-p 000b8000 08:06 656444     /usr/lib/i386-linux-gnu/libgfortran.so.3.0.0
00254000-00255000 rw-p 00000000 00:00 0 
0043a000-00594000 r-xp 00000000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00594000-00595000 ---p 0015a000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00595000-00597000 r--p 0015a000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00597000-00598000 rw-p 0015c000 08:06 658988     /lib/i386-linux-gnu/libc-2.13.so
00598000-0059b000 rw-p 00000000 00:00 0 
00714000-00738000 r-xp 00000000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00738000-00739000 r--p 00023000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00739000-0073a000 rw-p 00024000 08:06 659006     /lib/i386-linux-gnu/libm-2.13.so
00916000-00917000 r-xp 00000000 00:00 0          [vdso]
009fe000-00a18000 r-xp 00000000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00a18000-00a19000 r--p 00019000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00a19000-00a1a000 rw-p 0001a000 08:06 682158     /lib/i386-linux-gnu/libgcc_s.so.1
00eee000-00f0a000 r-xp 00000000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00f0a000-00f0b000 r--p 0001b000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
00f0b000-00f0c000 rw-p 0001c000 08:06 659007     /lib/i386-linux-gnu/ld-2.13.so
08048000-0804a000 r-xp 00000000 08:06 523755     /home/fluidistic/numerical2011/guide5/a.out
0804a000-0804b000 r--p 00001000 08:06 523755     /home/fluidistic/numerical2011/guide5/a.out
0804b000-0804c000 rw-p 00002000 08:06 523755     /home/fluidistic/numerical2011/guide5/a.out
09376000-09397000 rw-p 00000000 00:00 0          [heap]
b7600000-b7621000 rw-p 00000000 00:00 0 
b7621000-b7700000 ---p 00000000 00:00 0 
b7784000-b7786000 rw-p 00000000 00:00 0 
b77a3000-b77a5000 rw-p 00000000 00:00 0 
bff39000-bff5a000 rw-p 00000000 00:00 0          [stack]
Aborted
fluidistic@fluidistic-desktop:~/numerical2011/guide5$

When I open my output file (polynomial.dat), I see

   1.0000000      -0.0000000    
   2.0000000      -1.3000000    
   3.0000000      -0.0000000    
   4.0000000      -11.700000    
   5.0000000       0.0000000    
   6.0000000       487.50000

 

[Mark44]

Let's try some different data in the input file:
0 0
1 1
2 4

These are points on the graph of y = x^2

Try a step size of 0.5 and see how close you get to the correct values. That would give us an idea about the correctness of your implementation.

As far as the memory error, I don't think it is really an error. My guess it comes from you allocating memory for your arrays, but not deallocating that memory. Here's a link to a site that talks about array allocation and deallocation - http://www.stanford.edu/class/me200c/tutorial_90/07_arrays.html.

 

[fluidistic]

Let's try some different data in the input file:
0 0
1 1
2 4

These are points on the graph of y = x^2

Try a step size of 0.5 and see how close you get to the correct values. That would give us an idea about the correctness of your implementation.

As far as the memory error, I don't think it is really an error. My guess it comes from you allocating memory for your arrays, but not deallocating that memory. Here's a link to a site that talks about array allocation and deallocation - http://www.stanford.edu/class/me200c/tutorial_90/07_arrays.html.

Thanks for your concern with my program.
I modified kincaid.dat into your values, I chose n=2 and h=0.5.
I get as output file:

   0.0000000      -0.0000000    
  0.50000000     -0.12500000    
   1.0000000       0.0000000    
   1.5000000       1.1250000    
   2.0000000       4.0000000

Hmm...
Edit: Just tried with

1 1
2 4
3 9

as input file.
I got

1.0000000      -0.0000000    
   1.5000000     -0.18750000    
   2.0000000       0.0000000    
   2.5000000       1.6875000    
   3.0000000       6.0000000

as output.
The first p is always 0 and the third too.
Another Edit: Ah... this is because in the loop

do j=0, i-1
       prod=prod*(x_m-x(j))
    end do

, the first x_m-x(j), namely x_m-x(0) is worth 0. Hence when I multiply this by the next values for prod, I always get 0.
Hmm... how to fix this? I'm thinking.

 

[fluidistic]

News: I've got rid of the memory problem! I was assigning d an array from 1 to n. I modified it to 0 to n.
I realize that a loop is badly made. I start a do i=1, n instead of starting from 0.
But something in the algorithm provided doesn't make sense to me.
According to the book, $p(x)= \sum _{i=0}^n \prod _{j=0}^{i-1} (x-x_j)$. The first term doesn't make sense to me. If i=0, then j= what? j goes from 0 to... -1?!
I know the polynomial must look like $p(x)=d_0+d_1(x-x_0)+d_2(x-x_0)(x-x_1)+...$. So this is this part in the program I must modify, but I have no idea how.

 

[Mark44]

News: I've got rid of the memory problem! I was assigning d an array from 1 to n. I modified it to 0 to n.
I realize that a loop is badly made. I start a do i=1, n instead of starting from 0.
But something in the algorithm provided doesn't make sense to me.
According to the book, $p(x)= \sum _{i=0}^n \prod _{j=0}^{i-1} (x-x_j)$. The first term doesn't make sense to me. If i=0, then j= what? j goes from 0 to... -1?!
I know the polynomial must look like $p(x)=d_0+d_1(x-x_0)+d_2(x-x_0)(x-x_1)+...$. So this is this part in the program I must modify, but I have no idea how.

Here's your algorithm (with added indents).

for i=0, 1,...,n do
   di <- f(xi)
end
for j=1,2,...,n do
   for i=j, j+1, ..., n do
      di <-(di-d_{i-1} )/(xi-x_{i-j})
   end 
end.

And here's your code. It doesn't match very closely. You have a lot of extra loops that aren't in the algorithm. For example, you have a do while loop that isn't in the algorithm, and your inner do loops don't run from i = j to i = n.

If this were my code, I would rewrite it to follow the algorithm above.

subroutine Newton(x_m,h,x,d,y)
implicit none
integer:: i, j
real, intent(inout)::x_m,h
real:: prod, p
real, dimension(0:n), intent(inout) :: x,y
real, dimension(1:n), intent(inout) :: d
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
 do i=0,n
    d(i)=y(i)
! end do
  do j=1, n
     if (i>= j) then
    d(i)=(d(i)-d(i-1))/(x(i)-x(i-j))
     end if
  end do
end do
  prod=1.!x_m-x(0)
  do i=1, n  
    do j=0, i-1
       prod=prod*(x_m-x(j))
    end do
    p=d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine

 

[fluidistic]

Here's your algorithm (with added indents).

for i=0, 1,...,n do
   di <- f(xi)
end
for j=1,2,...,n do
   for i=j, j+1, ..., n do
      di <-(di-d_{i-1} )/(xi-x_{i-j})
   end 
end.

And here's your code. It doesn't match very closely. You have a lot of extra loops that aren't in the algorithm. For example, you have a do while loop that isn't in the algorithm, and your inner do loops don't run from i = j to i = n.

If this were my code, I would rewrite it to follow the algorithm above.

You are right. Even more than that, the given algorithm isn't even totally right. Because for each j, the d(i) are wrongly calculated. Not the first one but the second one is indeed wrong.
Because of the line "di <-(di-d_{i-1} )/(xi-x_{i-j})", d(i) uses d(i) to be redefined. In each j loop, d(i) changes more than once. Hmm I don't know how to explain. Instead of
di <-(di-d_{i-1} )/(xi-x_{i-j}), it should be di <-(yi-y_{i-1} )/(xi-x_{i-j}) which is totally different!
I've tested it with a friend (we lost 2 hours to spot this problem). Now we get a good d(i) generation. In other words, we do get right values for the coefficients d(i) of the interpolating polynomial.
I still have the problem of the $\prod$ loop. No matter what, prod=0 due to the definition of p(x) given in the book.
The definition doesn't make sense to me. $p(x)= \sum _{i=0}^n \prod _{j=0}^{i-1} (x-x_j)$. So for i=0, j=0 I get an x(-1). What does this even mean? x(-1) doesn't even exist, how can Fortran deal with this?
Also in the subroutine I had forgotten to define n.
Here's my current code for it:

subroutine Newton(x_m,h,x,d,y,n)

implicit none
integer,intent(in)			:: n
integer					:: i, j
real, intent(inout)			:: x_m,h
real					:: prod, p
real, dimension(0:n), intent(inout)	:: x,y,d

do i=0,n
d(i)=y(i)
end do
do j=1,n
  do i=j,n
  d(i)=(y(i)-y(i-1))/(x(i)-x(i-j))
  end do
    do i=j,n
    y(i)=d(i)
    end do 
end do
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
  prod=1.
  do i=0,n  
    do j=0,i-1
    prod=prod*(x_m-x(j))
    end do
    p=d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine
end program

Feel free to give me an insight on the definition of p(x) which doesn't make any sense to me.

 

[Mark44]

You are right. Even more than that, the given algorithm isn't even totally right. Because for each j, the d(i) are wrongly calculated. Not the first one but the second one is indeed wrong.
Because of the line "di <-(di-d_{i-1} )/(xi-x_{i-j})", d(i) uses d(i) to be redefined.

I suspect that this is by design. Each time d(i) is recalculated, based on a previous value of d(i). I could be wrong.

In each j loop, d(i) changes more than once. Hmm I don't know how to explain. Instead of
di <-(di-d_{i-1} )/(xi-x_{i-j}), it should be di <-(yi-y_{i-1} )/(xi-x_{i-j}) which is totally different!
I've tested it with a friend (we lost 2 hours to spot this problem). Now we get a good d(i) generation. In other words, we do get right values for the coefficients d(i) of the interpolating polynomial.
I still have the problem of the $\prod$ loop. No matter what, prod=0 due to the definition of p(x) given in the book.
The definition doesn't make sense to me. $p(x)= \sum _{i=0}^n \prod _{j=0}^{i-1} (x-x_j)$. So for i=0, j=0 I get an x(-1). What does this even mean? x(-1) doesn't even exist, how can Fortran deal with this?

The product notation makes sense only when i - 1 (the upper limit) is >= 0; i.e., when i is >= 1. When i = 1, you get one factor: (x - x₀); when i = 2, you get two factors: (x - x₀)(x - x₁). And so on.

Also in the subroutine I had forgotten to define n.

I noticed this before. The subroutine was using the value of n from the main program. IOW, n was a global variable. Passing n as a parameter is a better idea.

I'll try to take a look at your latest code a little later.

 

[fluidistic]

I suspect that this is by design. Each time d(i) is recalculated, based on a previous value of d(i). I could be wrong.

I find it hard to explain as I myself probably doesn't understand well. We did some algebra on a paper to see what the given algorithm did and it was wrong. For a given j, i will go from 0 to n and calculate several d(i)'s. The problem is that these d(i) must refer to old d(i). Not the ones the algorithm calculate for each j but the ones assigned with f(x_i).
We've tested on 2 polynomials (degree 3 and 4) in which we knew what were worth the coefficients d(i)'s and the program gave good results only with my/our new algorithm where we refer to y(i)'s.

The product notation makes sense only when i - 1 (the upper limit) is >= 0; i.e., when i is >= 1. When i = 1, you get one factor: (x - x₀); when i = 2, you get two factors: (x - x₀)(x - x₁). And so on.

Ah I see! Thanks a bunch for the info. I actually tried to fix it in my program but without success. I'm out for today (past 1 am) but will come back to it tomorrow for sure.

I noticed this before. The subroutine was using the value of n from the main program. IOW, n was a global variable. Passing n as a parameter is a better idea.

I'll try to take a look at your latest code a little later.

About n: I realized this since for a weird reason I got an error when compiling (I had no error before...) so I fixed it.

Thanks again for your concern, I appreciate it very much. Have fun :)

 

[fluidistic]

Now it works great! My friend helped me for it.
I had an error with the p=... line.
I modified a bit the look of my program (I like useless huge number for the character statement lol. I also like very condensed program, hardly readable for a human eye), here it is if you're (or anyone else) interested.

program div
implicit none
character(1000000)::datos
real(8), dimension(:), allocatable:: x,y,d
integer::n,i
real(8)::h,x_m

write(*,*)"Enter the degree n of the polynomial and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file containing the data to be interpolated"
read(*,'(A)')datos
allocate(x(0:n),y(0:n), d(0:n))
open(1,file=datos)
do i=0, n
read(1,*)x(i),y(i)
end do
close(1)
x_m=x(0)
call Newton(x_m,h,x,d,y,n)
contains
subroutine Newton(x_m,h,x,d,y,n)
implicit none
integer,intent(in)::n
integer::i,j
real(8), intent(inout)::x_m,h
real(8)::prod,p
real(8), dimension(0:n), intent(inout)::x,y,d

do i=0,n
d(i)=y(i)
end do
do j=1,n
do i=j,n
d(i)=(y(i)-y(i-1))/(x(i)-x(i-j))
end do
do i=j,n
y(i)=d(i)
end do 
end do
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
p=d(0)
  do i=1,n
    prod=1d0
    do j=0,i-1
    prod=prod*(x_m-x(j))
    end do
    p=p+d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine
end program

It interpolates very well the quadratic function. Notice that I've used double precision (8) for this final version.
Problem solved! Thank you for all your help and time.

 

[Mark44]

Now it works great! My friend helped me for it.
I had an error with the p=... line.
I modified a bit the look of my program (I like useless huge number for the character statement lol.

This could come back and bite you by causing your program to run out of memory. You are using a million bytes to store the name of the input file. It's a good idea to allocate a little extra storage, of course, but a million bytes of wasted storage isn't a good idea.

I also like very condensed program, hardly readable for a human eye)

I hope you're kidding. Some programmers take pride in writing dense, almost unreadable code, but in the real world, much of the life of a program is in maintenance, where a different programmer has to spend time fixing bugs that the original programmer introduced by being overly clever.

If program code is difficult to read, it will be even more difficult to understand, which makes it much more difficult to fix or change.

, here it is if you're (or anyone else) interested.

One thing you can do to make your code more readable is to indent the bodies of loops. I've done this in your code below.

program div
implicit none
character(1000000)::datos
real(8), dimension(:), allocatable:: x,y,d
integer::n,i
real(8)::h,x_m

write(*,*)"Enter the degree n of the polynomial and a step-size h" 
read(*,*)n,h
write(*,*)"Enter the name of the file containing the data to be interpolated"
read(*,'(A)')datos
allocate(x(0:n),y(0:n), d(0:n))
open(1,file=datos)
do i=0, n
   read(1,*)x(i),y(i)
end do
close(1)
x_m=x(0)

call Newton(x_m,h,x,d,y,n)
contains
subroutine Newton(x_m,h,x,d,y,n)
implicit none
integer,intent(in)::n
integer::i,j
real(8), intent(inout)::x_m,h
real(8)::prod,p
real(8), dimension(0:n), intent(inout)::x,y,d

do i=0,n
   d(i)=y(i)
end do
do j=1,n
   do i=j,n
      d(i)=(y(i)-y(i-1))/(x(i)-x(i-j))
   end do

   do i=j,n
      y(i)=d(i)
   end do 
end do
open(unit=2, file="polynomial.dat", status='old')
do while (x_m <= x(n))
   p=d(0)
   do i=1,n
      prod=1d0
      do j=0,i-1
         prod=prod*(x_m-x(j))
      end do
      p=p+d(i)*prod
  end do    
  write(2,*)x_m, p
  x_m=x_m+h
end do
close(2)
end subroutine
end program

It interpolates very well the quadratic function. Notice that I've used double precision (8) for this final version.
Problem solved! Thank you for all your help and time.

 

[fluidistic]

This could come back and bite you by causing your program to run out of memory. You are using a million bytes to store the name of the input file. It's a good idea to allocate a little extra storage, of course, but a million bytes of wasted storage isn't a good idea.

Oh I didn't know. In this case I'm going to choose 20 instead of 1 million. (I'm still laughing about my ignorance!)

I hope you're kidding. Some programmers take pride in writing dense, almost unreadable code, but in the real world, much of the life of a program is in maintenance, where a different programmer has to spend time fixing bugs that the original programmer introduced by being overly clever.

If program code is difficult to read, it will be even more difficult to understand, which makes it much more difficult to fix or change.

I'm not that kidding. But once again I think it would be indeed better if it was more readable. I have the feeling that if a program is dense, then it will be executed faster. Is it true? Or maybe it will compile faster but execute at the same speed than the same program but less dense.
Say I write

"program helloworld
(10^9 blank lines+ useless comments)
write(*,*)"Hey there"
end program"

Would the program compile as fast as if there was no blank lines+comments? What about its execution's speed?

One thing you can do to make your code more readable is to indent the bodies of loops. I've done this in your code below.

Thanks once again, I'm definitely going to use it this way.

 

[SteamKing]

Compilers strip out comments and white space (blank spaces/empty lines) before starting the lexical analysis of the program text. The incremental time spent by the compiler doing this is small compared to the time a programmer would spend trying to make sense of densely coded program code.

Besides, the time a good program takes to compile and link should allow the programmer to get a cup of coffee or take a stretch.

 

[fluidistic]

Compilers strip out comments and white space (blank spaces/empty lines) before starting the lexical analysis of the program text. The incremental time spent by the compiler doing this is small compared to the time a programmer would spend trying to make sense of densely coded program code.

Besides, the time a good program takes to compile and link should allow the programmer to get a cup of coffee or take a stretch.

Thanks for the information. Actually I've found almost by error (I was looking for Fortran codes on google) and found http://www.nikhef.nl/~templon/fortran/fortran_style.
Basically the guy states some " Fortran Coding Style for the Modern Programmer".
Number 5 is

Another good reason to leave out comments is that they cause slower
execution of your program.

And the golden rule is

17. The Golden Rule:
Regardless of the capacity of the machine your code is to run on,
it is more important to make it small and run blindingly fast,
than to maximize program readability and maintainability.

So... now I'm really confused. It seems like I was "right" thinking that a very dense/compact hard for the human eye to read code is noticeably "better" for the computer??

 

[SteamKing]

The guy who posted "Fortran Coding Style for the Modern Programmer" doesn't know how compilers work. Comments are put into a program file for the benefit of those who must understand and maintain the code (often not the original author/developer). Comments by nature don't mean anything to the computer executing the program and are left out of the code files generated by the compiler. Similarly, white spaces (blank spaces/lines) used to add clarity to a program text or illustrate program structure also don't mean anything to the computer executing the program. These too are omitted.

The speed of a compiled program depends on many factors, including these: the speed of the processor, whether the most efficient algorithm is used by the programmer, the amount of memory available, and how well the compiler generates efficient, optimized code fro the processor to execute. The compiler can't control how fast the processor is or how much memory is available, nor can it guarantee that a programmer writes only the best program. The compiler can only generate code based on the variables and commands contained in the programs it compiles. Different compilers produce code for the same program text which may be larger or smaller in size, or which may execute slower or faster on the same computer. Since the compiler itself is a sophisticated computer program, its efficiency and the quality of the code it produces often depends on the skill of the programmer/programmers who wrote it in the first place.

 

[Mark44]

Thanks for the information. Actually I've found almost by error (I was looking for Fortran codes on google) and found http://www.nikhef.nl/~templon/fortran/fortran_style.
Basically the guy states some " Fortran Coding Style for the Modern Programmer".
Number 5 is

Another good reason to leave out comments is that they cause slower
execution of your program.

And the golden rule is

17. The Golden Rule:
Regardless of the capacity of the machine your code is to run on,
it is more important to make it small and run blindingly fast,
than to maximize program readability and maintainability.

So... now I'm really confused. It seems like I was "right" thinking that a very dense/compact hard for the human eye to read code is noticeably "better" for the computer??

The guy who posted "Fortran Coding Style for the Modern Programmer" doesn't know how compilers work.

I'm pretty sure Glen Reesor actually knows what he is talking about, but was writing tongue-in-cheek. It's not very clear that he's pulling our leg, except that every item in his list is something you SHOULD NOT do. The reasons he gives for writing the list give us a clue about his actual intentions.

In my (relatively) short time modifying and working with legacy Fortran
codes, I have come across a number of, shall we say, interesting tendencies.
I attribute these tendencies to a number of things:

1. Historical limitations of compilers and hardware (very understandable--
at least for dusty-deck code, but depressing for newly written code).

2. The author being someone with zero training in how to write
maintainable, readable, quality software. (In my opinion, a very
significant factor in new Fortran software being written.)

Although, I keep telling myself that item 2 is very common, I've
lost count of the number of times I've been trying to decipher
a code segment and decided there should be an item 3:

3. The author of the code is from a different planet where their logic
is a complete reversal (yet twisted in some way) to ours :-)

 

[SteamKing]

Is Glen Reesor Scandinavian by any chance? There's another thread in this forum about Scandinavian programmers.