Function problem

Tags:
1. Aug 31, 2015

Galizius

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. Aug 31, 2015

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.

3. Aug 31, 2015

BvU

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

4. Aug 31, 2015

Staff: Mentor

I just noticed another problem:

The code is accessing x(0) and y(0), which are outside the array boundaries.

5. Aug 31, 2015

Staff: Mentor

I guess you didn't mean dividing by 0, in which case you beat me to it!

6. Aug 31, 2015

Galizius

That solved my problem! Thank You so much, I do not know how could I miss that...

7. Aug 31, 2015

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.

8. Aug 31, 2015

Galizius

Yes, I already did. Thank You again!

9. Aug 31, 2015