forums.silverfrost.com Forum Index forums.silverfrost.com
Welcome to the Silverfrost forums
 
 FAQFAQ   SearchSearch   MemberlistMemberlist   UsergroupsUsergroups   RegisterRegister 
 ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

subtle error in code?

 
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General
View previous topic :: View next topic  
Author Message
ognik



Joined: 04 Feb 2015
Posts: 19

PostPosted: Sun Feb 08, 2015 1:34 am    Post subject: subtle error in code? Reply with quote

As part of my course I have written program to evaluate a complex integral numerically using Simpsons rule. To avoid singularities from the limits, I split the integral into parts, each with its own equivalent limits. I coded each of these parts as a function. I believe my maths is right, but the program is not converging as it should - I have spent days checking both. I suspect that due to my newness to fortran, I have a subtle error in my code and would appreciate some help finding it, I've reached the stage I have looked at it too often and a fresh pair of more experienced eyes would be a good idea.
(sorry the code is not indented, I don't know how to preserve the indentation when copying)
---------------------------------
program simpson
!Summary: Integral from a to b of f(x) = h/3( f(a) + 4*f(a+h) + 2*f(a+2h) + 4*f(a+3h) + ...4*f(b-h) + f(b))
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
real, parameter :: pi = 3.1415927
real (kind=dp) :: h1, h2, part1, part2, exact, approx, diff, FUNC_1, FUNC_2, a, b, c, d, arg, temp
integer, parameter :: nn=50
integer :: N, cntr, coord
character(len=nn) :: stars=REPEAT('*',nn)

!----- Initialise ----
exact = (2*pi)/sqrt(3.0) ! Given
a=0.0**(1/3.0)
b=0.5**(1/3.0) ! 1st function integrates from 0 to 0.5 - change of variable u=t(1/3)
d=(1.0-0.5)**(3/2.0) ! 2nd function from 0.5 to 1.0 - after change of variable v=(1-t)**(3/2),
c=(1.0-1.0)**(3/2.0) ! - swap limits and change sign of integral
!----- header ----
print *, " "
print *, 'Program Simpson starts'
print *, 'Simpson approximation of Integral for f(t) = (t**(-2/3)) * ((1-t)**(-1/3))'
print *, 'Exact value = ', exact
print *, stars
print *, 'Enter value of N (multiples of 4) (0 to stop)'
print 20, 'N approximate ','difference '
read (*,*) N
do while ((mod(N,4).ne.0).OR.(n.lt.0))
print *, 'N must be positive and a multiple of 4 (or 0) - please re-enter N'
read *, N
end do
!----- end header ----
do while (N.ne.0) ! allows repeated trials of diff. values of h
! Evaluate 1st part
h1=(b-a)/N ! N intervals of h between a and b
print 10, ' Evalute first part from a= ', a, ' to b= ', b, ' at intervals h= ', h1 !debug
print *, " "
arg=a
approx=FUNC_1(arg) !1st term, f0
Do cntr = 1,N-1
if (mod(cntr,2)>0) then !set co-ord - either 4 (even) or 2 (odd)
coord = 2
else
coord = 4
end if
arg=arg+h1 !co-ord terms, f(a+h)
temp=FUNC_1(arg)
approx=approx+(coord*temp)
! print 40, ' term= ', cntr, ' coord= ', coord, ' arg= ', arg, ' function= ', temp, ' approx= ',approx !debug
end do
part1 = (h1/3.0)*(approx + FUNC_1(b)) !last term + simpson h/3 multiplier
! Evaluate 2nd part
h2=(d-c)/N ! N intervals of h between c and d
print 10, ' Evalute 2nd part from c= ', c, ' to d= ', d, ' at intervals h= ', h2
arg=c
approx=FUNC_2(arg) !1st term
Do cntr = 1,N-1
if (mod(cntr,2)>0) then !set co-ord - either 4 (even) or 2 (odd)
coord = 2
else
coord = 4
end if
arg=arg+h2
temp=FUNC_2(arg)
approx=approx+(coord*temp)
! print 40, ' term= ', cntr, ' coord= ', coord, ' arg= ', arg, ' function= ', temp, ' approx= ',approx !debug
end do
part2=(h2/3.0)*(approx + FUNC_2(d)) !last term + simpson h/3 multiplier
print *, part1, ' ', part2 !debug
approx = part1 + part2
diff=exact-approx
print 30, approx, diff
read (*,*) N
do while ((mod(N,4).ne.0).OR.(n.lt.0))
print *, 'N must be positive and a multiple of 4 (or 0) - please re-ent
Back to top
View user's profile Send private message
PaulLaidler
Site Admin


Joined: 21 Feb 2005
Posts: 7916
Location: Salford, UK

PostPosted: Sun Feb 08, 2015 10:01 am    Post subject: Reply with quote

Have you tried using the debugger SDBG and looking at the results on a line-by-line basis as they are evaluated. If necessary create intermediate values to confirm some calculations. Be careful with mixed integer/real expressions particularly when dividing. Use 1.0 for real 1 etc. Avoid dividing by small values.
Back to top
View user's profile Send private message AIM Address
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Sun Feb 08, 2015 11:35 am    Post subject: Reply with quote

The integral may be recognized as the Euler Beta function B(1/3, 2/3), see http://en.wikipedia.org/wiki/Beta_function.

Apart from losing the indentation, your program is also incomplete, so it is not possible to comment. Please repost, select the entire block of code with the mouse, and then click on the Code button. Or, if the code is too long, post it to the cloud with public visibility and provide a link to that location.
Back to top
View user's profile Send private message
ognik



Joined: 04 Feb 2015
Posts: 19

PostPosted: Sun Feb 22, 2015 2:37 am    Post subject: Reply with quote

Sorry for the delay, didn't get a notification about these replies.
I have pasted in the full code below...
[code:1:bcc4281167] program simpson
!Summary: Integral from a to b of f(x) = h/3( f(a) + 4*f(a+h) + 2*f(a+2h) + 4*f(a+3h) + ...4*f(b-h) + f(b))
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
real, parameter :: pi = 3.1415927
real (kind=dp) :: h1, h2, part1, part2, exact, approx, diff, FUNC_1, FUNC_2, a, b, c, d, arg, temp
integer, parameter :: nn=50
integer :: N, cntr, coord
character(len=nn) :: stars=REPEAT('*',nn)

!----- Initialise ----
exact = (2*pi)/sqrt(3.0) ! Given
a=0.0**(1/3.0)
b=0.5**(1/3.0) ! 1st function integrates from 0 to 0.5 - change of variable u=t(1/3)
d=(1.0-0.5)**(3/2.0) ! 2nd function from 0.5 to 1.0 - after change of variable v=(1-t)**(3/2),
c=(1.0-1.0)**(3/2.0) ! - swap limits and change sign of integral
!----- header ----
print *, " "
print *, 'Program Simpson starts'
print *, 'Simpson approximation of Integral for f(t) = (t**(-2/3)) * ((1-t)**(-1/3))'
print *, 'Exact value = ', exact
print *, stars
print *, 'Enter value of N (multiples of 4) (0 to stop)'
print 20, 'N approximate ','difference '
read (*,*) N
do while ((mod(N,4).ne.0).OR.(n.lt.0))
print *, 'N must be positive and a multiple of 4 (or 0) - please re-enter N'
read *, N
end do
!----- end header ----
do while (N.ne.0) ! allows repeated trials of diff. values of h
! Evaluate 1st part
h1=(b-a)/N ! N intervals of h between a and b
print 10, ' Evalute first part from a= ', a, ' to b= ', b, ' at intervals h= ', h1 !debug
print *, " "
arg=a
approx=FUNC_1(arg) !1st term, f0
Do cntr = 1,N-1
if (mod(cntr,2)>0) then !set co-ord - either 4 (even) or 2 (odd)
coord = 2
else
coord = 4
end if
arg=arg+h1 !co-ord terms, f(a+h)
temp=FUNC_1(arg)
approx=approx+(coord*temp)
! print 40, ' term= ', cntr, ' coord= ', coord, ' arg= ', arg, ' function= ', temp, ' approx= ',approx !debug
end do
part1 = (h1/3.0)*(approx + FUNC_1(b)) !last term + simpson h/3 multiplier
! Evaluate 2nd part
h2=(d-c)/N ! N intervals of h between c and d
print 10, ' Evalute 2nd part from c= ', c, ' to d= ', d, ' at intervals h= ', h2
arg=c
approx=FUNC_2(arg) !1st term
Do cntr = 1,N-1
if (mod(cntr,2)>0) then !set co-ord - either 4 (even) or 2 (odd)
coord = 2
else
coord = 4
end if
arg=arg+h2
temp=FUNC_2(arg)
approx=approx+(coord*temp)
! print 40, ' term= ', cntr, ' coord= ', coord, ' arg= ', arg, ' function= ', temp, ' approx= ',approx !debug
end do
part2=(h2/3.0)*(approx + FUNC_2(d)) !last term + simpson h/3 multiplier
print *, part1, ' ', part2 !debug
approx = part1 + part2
diff=exact-approx
print 30, approx, diff
read (*,*) N
do while ((mod(N,4).ne.0).OR.(n.lt.0))
print *, 'N must be positive and a multiple of 4 (or 0) - please re-enter N'
read *,
Back to top
View user's profile Send private message
DanRRight



Joined: 10 Mar 2008
Posts: 2813
Location: South Pole, Antarctica

PostPosted: Sun Feb 22, 2015 3:58 pm    Post subject: Reply with quote

Paul & Co, How about introducing temporal storage folder (say, keeping files for a year or so) on this site to accommodate larger source examples? The limitation like this forum has right now which was introduced by phpBB 10-15 years ago when even servers worked on 10 GB harddrives maximum looks very archaic. Cutting files on small pieces introduces errors and is an unjustified headache for almost everyone. What is the total size of all posts on this entire forum, does it really need too much support and huge drives for backup?
Back to top
View user's profile Send private message
silverfrost
Site Admin


Joined: 29 Nov 2006
Posts: 191
Location: Manchester

PostPosted: Mon Feb 23, 2015 11:49 am    Post subject: Reply with quote

I was hoping there would be an' enable attachments' check-box in the forum admin section but I cannot find one at the moment. An alternative is to use a gist:

https://gist.github.com/silverfrost/fe90e13e9d1b1054f2d3
Back to top
View user's profile Send private message Visit poster's website
Wilfried Linder



Joined: 14 Nov 2007
Posts: 314
Location: Düsseldorf, Germany

PostPosted: Mon Feb 23, 2015 6:35 pm    Post subject: Reply with quote

In the administrator's area, select the first tab ("Common" ?), then "Board configuration", then "file attachments", now you can enable attachments and set their maximum size. (I only know the German expressions in a phpBB forum, so the English ones may differ from what I wrote).

Wilfried
Back to top
View user's profile Send private message
silverfrost
Site Admin


Joined: 29 Nov 2006
Posts: 191
Location: Manchester

PostPosted: Mon Feb 23, 2015 6:56 pm    Post subject: Reply with quote

It isn't in there and a search through the source doesn't throw anything up either. It is quite an old versioj of phpBB, maybe it isn't supported.
Back to top
View user's profile Send private message Visit poster's website
ognik



Joined: 04 Feb 2015
Posts: 19

PostPosted: Sun Mar 01, 2015 10:55 pm    Post subject: Code posted Reply with quote

Hi - apologies, the code was complete in preview, chopped off with submit.

I have posted the program here - https://drive.google.com/file/d/0ByNoaXX7FNqRSVpfY1hKekR4Vk0/view?usp=sharing, click and it should display correctly in a browser.

I have run this with debugger and the program gets the manually calculated values (the use of extra variables was mostly for debugging) - but too different from the exact value I have for testing the algorithm. Gone over it all many times and just not getting anything new ....
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Mon Mar 02, 2015 2:44 am    Post subject: Reply with quote

The error is not in your programming but in your mathematical work. Specifically, check the exponents in the expressions for the integration limits c and d, and in the definition of v. Your intention probably is that, near t=1, the integrand is approximated by 1/(1-t)^(1/3), and you should define v such that dv = -dt/(1-t)^(1/3).

By presenting the program in a cluttered way, you made it uninviting for people to try to help you. Furthermore, instead of describing the integral in a straightforward mathematical way, you forced people to infer that from your code.

Communication, communication!
Back to top
View user's profile Send private message
ognik



Joined: 04 Feb 2015
Posts: 19

PostPosted: Wed Mar 18, 2015 10:32 pm    Post subject: Maths checked Reply with quote

Hi - sorry I took a while to re-check the math & program and get back (and I have a day job Smile). Still stuck unfortunately.....

The problem I'm trying to program is the integral wrt t, between 0 and 1, of f(t)=(t**(-2/3))*(1-t)**(-1/3).
Both the limits give singularities for this function, so I have to split the integral into 2 parts, each part deals with 1 of the singularities. I arbitrarily split them at 0 -> 1/2, 1/2 -> 1

Part A) Let u=t**(1/3). Then 3du=t**(-2/3)dt and t=u**3
Limits: 0 .LE. t .LE. 1/2, gives 0 .LT. u .LT. (1/2)**(1/3)

So 1st part becomes Integral wrt u, between (a in the program) 0 and (b) (1/2)** (1/3), of g(u)=3*(1-u**3)**(-1/3)

B) Let v=(1-t)**(2/3). Then -3/2dv=(1-t)**(-1/3)dt and t=1-v**(3/2)
Limits (1/2)**(3/2) .LE. t .LE.1, gives (1/2)**(3/2) .LE. v .LE. 0

So 2nd part becomes integral wrt v, of g(v) =(-3/2)*(1-v**(3/2))**(-2/3)
I then swapped the limites to remove the - in front of the integral, so limits are c=0 to d=(1/2)**(3/2)
------------
Then the progfram should estimate the 2 integrals separately, using simpsons method.
h is the step size, I use the same number of steps in each part, which means a different step size; I'm pretty sure this is OK.
Some of my coding is from me trying to find the error, so its become a bit verbose to help me check through it using debugger - the program seems to do exactly what I expect from manual calcs - so after many hours at this I just can't see where I have gone wrong - math or program.....
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2554
Location: Sydney

PostPosted: Wed Mar 18, 2015 10:52 pm    Post subject: Reply with quote

Quote:
The problem I'm trying to program is the integral wrt t, between 0 and 1, of f(t)=(t**(-2/3))*(1-t)**(-1/3).


You could use Newton's rule or Simpsons rule, however there is an alternative approach, which is simply getting an unbiased average of the function between 0 - 1.
This can be achieved using a random number generator and taking the average of the function results for this random average. You can first try:
Code:
s = 0
do I = 1,n
 call radom_number (x)
 s = s + f(x)
end do
integral = s / n

You could then modify to replace "n" with a convergence test.
As N gets larger, you may need to correct for round-off when calculating S + f(x). There are numerical methods for accumulating and applying the error at each iteration. ( I have forgotten the algorithm's name!)
The point I am making is that this is a valid alternative approach with a controlled accuracy.

John

The following is a program I wrote to test the accuracy estimate. I'm sure it can be improved.
Code:
real*8 function f (t)
  real*8 t
  real*8, parameter :: power = 1.0d0/3.0d0
!
  f = (t/(1-t))**power/t
end function f

subroutine do_integral (n_i)
  integer*4 n_i, n_m, i, m
  real*10   s
  real*8    x, integral
  real*8    sum_i, sum_ii, av, sd
  real*8    f
  external  f
!
  n_m = 10
!
  sum_i  = 0
  sum_ii = 0
  do m = 1,n_m
    s = 0
    do I = 1,n_i
      call random_number (x)
      s = s + f(x)
    end do
    integral = s / dble(n_i)
    sum_i    = sum_i  + integral
    sum_ii   = sum_ii + integral**2
    write (*,*) m, integral
  end do
!
  x  = n_m
  av = sum_i/x
  sd = sqrt ( sum_ii/x - av**2 )
  write (*,*) ' '
  write (*,*) 'number of integrals   =', n_m
  write (*,*) 'samples per integral  =', n_i
  write (*,*) 'average  of integrals =', av
  write (*,*) 'Stdev of integrals    =', sd
  write (*,*) 'Stdev of average      =', sd / sqrt ( x )
!
end subroutine do_integral

!  program to generate integral
!
  integer*4 n,i
  n = 1
  do i = 1,7
    n = n*10
    call do_integral (n)
  end do
end


Last edited by JohnCampbell on Thu Mar 19, 2015 1:02 am; edited 2 times in total
Back to top
View user's profile Send private message
mecej4



Joined: 31 Oct 2006
Posts: 1885

PostPosted: Wed Mar 18, 2015 11:13 pm    Post subject: Reply with quote

Ognik, you have repeated the mistake. Look at the line in your post of today that begins with "B)". Substitute t = 1/2 in the expression for v. Compare that with what you wrote as the upper limit of v. [spoiler]look at the exponent.[/spoiler]

Aside: I find it hard to believe that there are educational institutions today that give you over a month to work out a homework exercise. Some of us may remember that in the old days (before there were computers) we would be asked to do the calculation with the help of log tables, and do so in less than a week.

John, your suggested Monte Carlo approach can be quite instructive, but in this case, where the integrand is singular at both limits, the results can be expected to be poor when the original variables are used. The transformed integrals, on the other hand, will probably behave fine with Monte Carlo. However, I don't know if Ognik has been given free choice as to the algorithm that he may use. And, but for the tiny slip, he is almost there with a working program.
Back to top
View user's profile Send private message
JohnCampbell



Joined: 16 Feb 2006
Posts: 2554
Location: Sydney

PostPosted: Thu Mar 19, 2015 1:46 am    Post subject: Reply with quote

There is a method of reducing round-off error, reported by Kahan.
It can be coded as:
Code:
real*8 function f (t)
!
!  Calculates the function value inside 0:1
  real*8 t, tm
  real*8, parameter :: power = -1.0d0/3.0d0
!
! old  f = (t/(1-t))**(1.d0/3.d0) / t
!
  f  = 0
  tm = 1 - t
  if ( t > 0 .and. tm > 0 ) then
     f  = (t*t*tm)
     if ( f > 0 ) f = f ** power
  end if
end function f

subroutine do_integral (n_i)
  integer*4 n_i, n_m, i, m
  real*8    s, s1
  real*8    x, integral
  real*8    sum_i, sum_ii, av, sd
  real*8    f, fx, e, em, fxm
  real*8    pi, exact
  external  f
!
  pi = 4 * atan (1.0d0)
  exact = 2 * pi / sqrt (3.0d0)
!
  n_m = 15
!
  sum_i  = 0
  sum_ii = 0
  do m = 1,n_m
    s  = 0
    e  = 0
    em = 0
    fxm = 0
    do i = 1,n_i
      call random_number (x)
!      s = s + f(x)
!
!   Alternative to accumulate and return error (Kahan algorithm)
      fx = f(x)
      s1 = s + (fx+e)
      e  = (fx+e) - (s1-s)
      em = max ( em, abs(e) )
      fxm = max ( fxm, fx )
      s  = s1
    end do
!
    integral = s / dble(n_i)
    sum_i    = sum_i  + integral
    sum_ii   = sum_ii + integral**2
    write (*,11) m, integral, s,   e,     em,    fxm 
!                   integral, sum, error, max e, max fx
11  format (i5,f15.10,f20.8,3es13.3)
  end do
!
  x  = n_m
  av = sum_i/x
  sd = sqrt ( (sum_ii - x*av**2)/(x-1) )
!
  write (*,*) ' '
  write (*,*) 'number of integrals   =', n_m
  write (*,*) 'samples per integral  =', n_i
  write (*,*) 'average  of integrals =', av
  write (*,*) 'Stdev of integrals    =', sd
  write (*,*) 'Stdev of average      =', sd / sqrt ( x-2 )
  write (*,*) 'Actual average Error  =', av - exact
!
end subroutine do_integral

The alternative of using higher precision real*10 for the accumulator appears to be more effective.
It is interesting to compare the error estimate to the actual error.

John
Back to top
View user's profile Send private message
ognik



Joined: 04 Feb 2015
Posts: 19

PostPosted: Thu Mar 26, 2015 2:05 am    Post subject: Working thanks Reply with quote

Thanks for your patience mecej4, the limit error was so obvious when I found it, I couldn't believe I kept missing it. Program now well behaved - and I learned a few extra things along the way. Thanks also John, although mecej4 is right in that I have no choice as to the method to use. But I got a fair accuracy using a high number of steps (~96)
Back to top
View user's profile Send private message
Display posts from previous:   
Post new topic   Reply to topic    forums.silverfrost.com Forum Index -> General All times are GMT + 1 Hour
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group