forums.silverfrost.com

[ognik]

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

[PaulLaidler]

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.

[mecej4]

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.

[ognik]

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 *,

[DanRRight]

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?

[silverfrost]

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

[Wilfried Linder]

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

[silverfrost]

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.

[ognik]

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 ....

[mecej4]

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!

[ognik]

Hi - sorry I took a while to re-check the math & program and get back (and I have a day job ). 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.....

[JohnCampbell]

[mecej4]

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.

[JohnCampbell]

There is a method of reducing round-off error, reported by Kahan.
It can be coded as:

[ognik]

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)