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

Function problem

  1. Aug 31, 2015 #1
    Hello I am trying to exchange following code into two functions which I will need to use in further programs.
    The code is as follows:

    Code (Text):
    do i=0,k
      do j=0,k
      if (i/=j) then
        lj(i)=lj(i)*((xx-x(i))/(x(j)-x(i)))
      end if
      enddo
      l(i)=l(i)+y(i)*lj(i)
      write(2,*) x(i),y(i),l(i)
    enddo
    I would like to swap lj and l arrays for the functions. I have tried something like that:

    Code (Text):
    real(8) function lj(j,xx,x,k)
      real(8) :: x(1000)
      real(8) :: xx
      integer :: i,j,k

      lj=1.d0
    do i=0,k
      do j=0,k
      if (i/=j) then
        lj=lj*((xx-x(i))/(x(j)-x(i)))
      end if
      enddo
    enddo
    Code (Text):
    real(8) function l(xx,x,y,k)
    real(8) :: x(1000),y(1000), xx
    integer :: i,j,k
      do i=0,k
        l=l+y(i)*lj(i,xx,x,k)
      enddo
    And then I am trying to call it in the main program
    Code (Text):
    do i=0,k
      write(2,*) x(i),y(i),l(x(i),x,y,k)
    enddo
    But it does not seem to be working as before. What am I doing wrong?
     
    Last edited: Aug 31, 2015
  2. jcsd
  3. Aug 31, 2015 #2

    DrClaude

    User Avatar

    Staff: Mentor

    This is very vague. You'll get better help when you clearly describe why (you think) it is wrong.

    l is not initialized before being used in the loop.
     
  4. Aug 31, 2015 #3

    BvU

    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Not sure. In my compiler arrays are 1-based and dividing by 0 is disastrous. From your snippets I can come to something that runs but I have no idea if it's what you intend:
    Code (Text):

    real(8) function lj(j,xx,x,k)
      real(8) :: x(1000)
      real(8) :: xx
      integer :: i,j,k
     
      lj=1.d0
     do i=1,k
      do j=1,k
      if (x(i) /= x(j)) then
      if (i/=j) then
      lj=lj*((xx-x(i))/(x(j)-x(i)))
      end if
      end if
      enddo
     enddo
    end
     real(8) function l(xx,x,y,k)
     real(8) :: x(*),y(*), xx,  lj
     external lj
     integer :: i,j,k
      do i=1,k
      l=l+y(i)*lj(i,xx,x,k)
      enddo
     end
     
    Program test
    real*8 l,x(1000),y(1000)
    external l
    k=4
    do i=1,k
      write(2,*) x(i),y(i),l(x(i),x,y,k)
    pause
    enddo
    end
     
    with output

    Code (Text):
      0.000000000000000E+000  0.000000000000000E+000  6.790386531088871E-313
    suspicious ! We should probably initialize a few things
     
  5. Aug 31, 2015 #4

    DrClaude

    User Avatar

    Staff: Mentor

    I just noticed another problem:

    The code is accessing x(0) and y(0), which are outside the array boundaries.
     
  6. Aug 31, 2015 #5

    DrClaude

    User Avatar

    Staff: Mentor

    I guess you didn't mean dividing by 0, in which case you beat me to it!
     
  7. Aug 31, 2015 #6
    That solved my problem! Thank You so much, I do not know how could I miss that...
     
  8. Aug 31, 2015 #7

    DrClaude

    User Avatar

    Staff: Mentor

    Be sure to also implement the fix to the problem I mentioned in post #2. It might work for now, but you never know when it might fail.
     
  9. Aug 31, 2015 #8
    Yes, I already did. Thank You again!
     
  10. Aug 31, 2015 #9

    Mark44

    Staff: Mentor

    A few comments about your code:
    1. The names you have picked for most of your variables are terrible! Single letter names such as i and j are OK for loop control variables, but names like lj, x, xx, and k don't give any indication of what you're going to use them for. There is no reason not to choose a name that suggest what the variable represents. Many years ago there were limits on how long a name you could use for a variable, so having variable names like x, xx, and lj is completely unwarranted, and makes your code harder to read than it should be.
    2. Why are you passing j as a parameter of your function? This variable gets set to 0, 1, 2, ..., k. It is strictly a local loop control variable, so there is no point in passing it as a function parameter.
    3. There isn't a single comment in your code.
    4. There isn't a single whitespace character in your code. The compiler doesn't care, but it's much harder for humans to read this -- lj=lj*((xx-x(i))/(x(j)-x(i))) -- than this --
      lj = lj * ( (xx - x(i) ) / ( x(j) - x(i) ) )
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook