Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Need help with fortran 90 on an exercice

  1. Apr 15, 2008 #1

    fluidistic

    User Avatar
    Gold Member

    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 [tex]P(x)=a_0+a_1(x-x_0)+...+a_n(x-x_o)...(x-x_{n-1})[/tex].
    And then it must output the value of the interpolating polynomial evaluated m times between [tex]x_0[/tex] and [tex]x_n[/tex]. 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 [tex]x_0, x_n[/tex]. m divided the interval [tex]x_0, x_n[/tex] 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 do


    Contains
    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
     
    Last edited: Apr 15, 2008
  2. jcsd
  3. Apr 15, 2008 #2

    alphysicist

    User Avatar
    Homework Helper

    Hi fluidistic,

    I believe I see two errors in your program. You calculate [itex]a_2 \to a_n[/itex] 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,

    [tex]a_2 = \frac{y_2 -a_1(x_2 - x_0) - a_0}{(x_2-x_0)(x_2-x_1)}[/tex]

    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 [itex]a_{k-1}[/itex] in your formula for [itex]a_k[/itex].
    ---------------------------------------------

    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 [itex]a_{n-1}[/itex]. 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).
     
  4. Apr 16, 2008 #3

    fluidistic

    User Avatar
    Gold Member

    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!
     
  5. Apr 16, 2008 #4

    alphysicist

    User Avatar
    Homework Helper

    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.
     
  6. Apr 16, 2008 #5
    the newton form of the interpolating polynomial is more difficult to code than the lagrange. i would suggest using the lagrange if you can.
     
  7. Apr 17, 2008 #6

    fluidistic

    User Avatar
    Gold Member

    Hello ice109,
    In my anterior exercise I had to do it with Lagrange's interpolation and it worked. For me both are difficult!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?



Similar Discussions: Need help with fortran 90 on an exercice
  1. Help with Fortran 90 (Replies: 0)

Loading...