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

Tricky Intregral for numerical quadrature

  1. Jan 21, 2015 #1
    Hi - I have just started 'Computational Physics' by Koonin & Meredith, - through distance learning.
    Exercise 1.3 needs a program to evaluate an integral - I'm finding myself kinda rusty on integrals. The hint says - split range of integration into parts, use different change of variable in each part... but I've spent most of a day trying different things (I suspect a trig substitution), so would appreciate a strong hint here - the actual exercise is to write the program :) Equation is ∫(t-2/3)(1-t)-1/3dt between 0 and 1.

  2. jcsd
  3. Jan 21, 2015 #2


    User Avatar
    Science Advisor

    The problem is that [itex]t^{-2/3}[/itex] goes to infinity as x goes to 0 and [itex](1- t)^{-1/3}[/itex] goes to infinity as x goes to 1. The suggestion is that you separate the interval into "0 to a" and "a to 1" where "a" is some number between 0 and 1, 1/2, say. Then, set x be equal to some function of t such that the function goes to infinity as t goes to 0.
  4. Jan 21, 2015 #3
    Hi, yes I was that far, my 'brick wall' is that I cant find any suitable substitutions to make - from which I would be able to figure out the transformed limits. I tried Integrating by parts and various substitutions, all of which complicated rather than simplified.... could you suggest a suitable split or substitution please?

  5. Jan 24, 2015 #4


    User Avatar
    Homework Helper


    $$t=\sin^3(x) \text{ or } t=\cos^3(x)$$

    that will remove the singularity on both sides
    you may or may not then wish to use
  6. Jan 24, 2015 #5
    I had tried t=sin3(x), probably I am making a mistake so would appreciate your checking as far as I got:
    dt/dx = 3sin2(x)cos(x) ⇒ dt=3sin2(x)cos(x)dx
    integral becomes ∫3cos(x) dx/(1-sin3(x))1/3)
    ...which still leaves me stumped....
  7. Jan 24, 2015 #6


    User Avatar
    Homework Helper

    That is right why are you stumped? The integral is very well behaved.
  8. Jan 24, 2015 #7
    Split integral into two, $$\int_0^1 t^{-2/3} (1-t)^{-1/3} dt = \int_0^{1/2} t^{-2/3} (1-t)^{-1/3} dt + \int_{1/2}^1 t^{-2/3} (1-t)^{-1/3} dt.$$ In the first integral the term ##(1-t)^{-1/3}## is bounded, so you only need to take care of the term ##t^{-2/3}##. Integrating ##\int t^{-2/3} dt = 3 t^{1/3} +C## you can see that using substitution ##u=t^{1/3}## gives you an integral of a bounded function and over a bounded interval, which you can then evaluate numerically. In fact, any substitution ##u = t^a##, ##0<a\le 1/3## works. but I feel the substitution ##u=t^{1/3}## works best.

    For the second integral you can act similarly, a substitution ##u =(1-t)^a##, ##0<a\le 2/3## works. Again I think ##a=2/3## would be best.
  9. Jan 25, 2015 #8
    Hi Lurflurf, it looks good to me too - until I try and resolve it. If it was 1-sin2 in the bracket or some simpler power ... then I could see a way forward, but I'm just fiddling without any intuition ... I will try Hankeye's approach tomorrow, getting late here - but thanks both.
  10. Jan 25, 2015 #9


    User Avatar
    Homework Helper

    as I mentioned above you can use

    $$\int_0^{\pi/2} \! \dfrac{\cos(x)}{\sqrt[3]{1-\sin^3(x)}}\, \mathrm{d}x=\int_0^{\pi/2} \! \sqrt[3]{\dfrac{\cos^3(x)}{1-\sin^3(x)}}\, \mathrm{d}x=\int_0^{\pi/2} \! \sqrt[3]{\cos(x)\dfrac{1+\sin(x)}{1+\sin(x)+\sin^2(x)}}\, \mathrm{d}x$$

    I am not sure that is better though
  11. Jan 25, 2015 #10


    User Avatar
    Homework Helper

    one other idea

    $$\int_0^1 \! t^{-2/3}(1-t)^{-1/3} \, \mathrm{d}t=\int_0^1 \! (t^{-2/3}+(1-t)^{-1/3}) \, \mathrm{d}t-\int_0^1 \! (t^{-2/3}+(1-t)^{-1/3}-t^{-2/3}(1-t)^{-1/3}) \, \mathrm{d}t$$
  12. Jan 25, 2015 #11
    Hi - Hawkeye's was easiest for me, please check my outcome:
    1st part: 3 ∫ (1 - u3))-1/3 du from u=0 to u= 0.125
    2nd part: -3/2 ∫ (1 - v3/2)-2/3 dv from v= 1 - (1/2)3/2 to 0.0
    I have checked this and it seems right, but would just appreciate confirmation - the program I wrote using this isn't giving the expected result...
  13. Jan 25, 2015 #12
    Just correcting my typos above, sorry:
    1st part: 3 ∫ (1 - u3)-/13 du - from u=0 to u= 0.125
    2nd part: -3/2 ∫ (1 - v3/2)-2/3 dv - from v= (1/2)3/2 to v = 0.0
    I have checked this and it seems right, but would just appreciate confirmation - the program I wrote using this isn't giving the expected result...
  14. Jan 26, 2015 #13


    User Avatar
    Homework Helper

    I got
    $$3\int_0^{2^{-1/3}} \! (1-u^{3})^{-1/3} \, \mathrm{d}u=\dfrac{\pi}{\sqrt{3}}+\log(2) \\ \dfrac{3}{2}\int_0^{2^{-2/3}} \! (1-u^{3/2})^{-2/3} \, \mathrm{d}u=\dfrac{\pi}{\sqrt{3}}-\log(2)$$
  15. Jan 26, 2015 #14
    Upper bound in the first part should be ##2^{-1/3}##, probably that is why you are getting an unexpected answer. The integrals should be as in the https://www.physicsforums.com/threads/tricky-intregral-for-numerical-quadrature.793526/members/lurflurf.19666/ [Broken]'s post above: try to use them in your program.

    PS. If you plug the integrals in Mathematica (or something similar) you get the answers as in the https://www.physicsforums.com/threads/tricky-intregral-for-numerical-quadrature.793526/members/lurflurf.19666/ [Broken]'s post. Or, if you plug the original integral, you get ##\frac{2\pi}{\sqrt 3}##. But the exact integration is quite complicated, the antiderivative of the original function involves hypergeometric functions, and I will not discuss it here.

    I understand that you do not need an exact answer, you need to compute the integral numerically. Applying standard numerical methods to integrals in the https://www.physicsforums.com/threads/tricky-intregral-for-numerical-quadrature.793526/members/lurflurf.19666/ [Broken]'s post you should get the result. The reason that you need to transform the integral before using a numerical method is that standard numerical methods do not work well when functions have singularities (when functions are unbounded).
    Last edited by a moderator: May 7, 2017
  16. Jan 26, 2015 #15
    Hi - I think I have the same bounds, first part with u=t1/3 I have the upper as (1/2)1/3 - which is the same as 2-1/3 - but the latter does look 'better'. Similarly the 2nd upper bound I have 1/23/2 which is the same as 2-3/2. I assume this will make no difference? ..... so - now to recheck the program....
  17. Jan 26, 2015 #16
    While I am checking the body - I am using the composite simpsons rule - here are the functions I coded for each part ...
    Function FUNC_1(u) !calc 3* [(1-(u**3))**(-1/3)]
    integer, parameter :: dp = selected_real_kind(15, 307)
    real(kind=dp) :: u, FUNC_1, temp, out
    ! print 10, 'u= ',u, ' temp= ', temp, ' out= ', out !debug
    Func_1= out
    10 format (3(A,F10.7))
    End Function

    Function FUNC_2(v) !calc (3/2)* [(1-(v**(3/2))**(-2/3)]
    integer, parameter :: dp = selected_real_kind(15, 307)
    real(kind=dp) :: v, FUNC_2, out
    ! print 10, 'v= ',v, ' out= ', out !debug
    Func_2= out
    10 format (2(A,F10.7))
    End Function
  18. Jan 30, 2015 #17
    Hi again, haven't managed to find the problem, my program looks OK but doesn't output values close enough to the actual answer to be correct; must be something wrong. Also it doesn't converge with smaller step size. The above functions are coded in fortran95, can anyone see an error in the way I have implemented the 2 integrals?
  19. Mar 18, 2015 #18
    Been going over both the math and the program and still can't find the error. To summarize the math:

    The problem I'm trying to solve numerically, is the integral wrt t, between 0 and 1, of f(t)=(t(-2/3))*(1-t)(-1/3).
    Both the limits give singularities, so I have to split the integral into 2 parts, each part deals with 1 of the singularities. 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=u3
    Limits: 0 ≤ t ≤ 1/2, gives 0 ≤ u ≤ (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-u3)(-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) ≤ t ≤1, gives (1/2)(3/2) ≤ v ≤ 0

    So 2nd part becomes integral wrt v, of g(v) =(-3/2)*(1-v(3/2))(-2/3)
    I then swapped the limits to remove the - in front of the integral, so limits are c=0 to d=(1/2)(3/2)
    Then the program should estimate the 2 integrals 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.....

    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.
  20. Mar 18, 2015 #19
    You got the upper limit in the second integral wrong, it should be ##(1/2)^{2/3}##, not ##(1/2)^{3/2}##, see the post # 13 for the correct integrals. And there are coefficient in front of each integral, did you include them? I could not see in your program where you add 2 integrals.
  21. Mar 18, 2015 #20
    Many thanks about the limit, dont know how I kept mixing that up. (I included the coefficients in the 2 functions (Eg out=3*out) )
    I've just tried the program and it behaves much better, never completely diverges though - the best error I get is 0.0000000026 for N=96 (once it was working I improved the exact calc with a full _DP value for pi).
    Do you know if the simpson method ever can converge exactly? (I suspect not)
  22. Mar 18, 2015 #21
    Usually not. But id your function is a polynomial of degree at most 3 (or even if it is a spline of degree at most 3, i.e.. on each interval ##[2k h, 2(k+1)h]## the function is polynomial of degree at most 3, different polynomials for different intervals) then the Sympson's method gives the exact answer. Otherwise, look for the estimate of the error here.
  23. Mar 19, 2015 #22
    Thanks again, you've been most helpful and - importantly - I clearly understand it all :-) .
    In the meantime I have got through ch 1 & 2 and bust with ch3 of Koonin; the sense I get is that there is a very broad set of algorithms for any particular numeric approach, each with some pro's and cons'. I wonder - is there a sort of cheat sheet somewhere that summarizes the 'best' few for each function type (derivative, integral, roots, ODE, PDE etc.)? The sort of parameters I am thinking about are error, complexity of programming, some sort of run-time index etc......
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Similar Threads for Tricky Intregral numerical
I Gaussian Quadrature on a Repeated Integral
I Derivative and Numerical approach