Tricky Intregral for numerical quadrature

[ognik]

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.

Thanks
ALan

 

[HallsofIvy]

The problem is that $t^{-2/3}$ goes to infinity as x goes to 0 and $(1- t)^{-1/3}$ 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.

 

[ognik]

Hi, yes I was that far, my 'brick wall' is that I can't 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?

Thanks
Alan

 

[lurflurf]

substitute

$$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
$$\frac{1-\sin^2(x)}{1-\sin^3(x)}=\frac{1+\sin(x)}{1+\sin(x)+\sin^2(x)}$$

 

[ognik]

I had tried t=sin³(x), probably I am making a mistake so would appreciate your checking as far as I got:
t=sin³(x)
dt/dx = 3sin²(x)cos(x) ⇒ dt=3sin²(x)cos(x)dx
integral becomes ∫3cos(x) dx/(1-sin³(x))^1/3)
...which still leaves me stumped...

 

[lurflurf]

That is right why are you stumped? The integral is very well behaved.

 

[Hawkeye18]

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.

 

[ognik]

Hi Lurflurf, it looks good to me too - until I try and resolve it. If it was 1-sin² 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.

 

[lurflurf]

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

 

[lurflurf]

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

 

[ognik]

Hi - Hawkeye's was easiest for me, please check my outcome:
1st part: 3 ∫ (1 - u³)^)-1/3 du from u=0 to u= 0.125
2nd part: -3/2 ∫ (1 - v^3/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...

 

[ognik]

Just correcting my typos above, sorry:
1st part: 3 ∫ (1 - u³)^-/13 du - from u=0 to u= 0.125
2nd part: -3/2 ∫ (1 - v^3/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...

 

[lurflurf]

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)$$

 

[Hawkeye18]

Just correcting my typos above, sorry:
1st part: 3 ∫ (1 - u³)^-/13 du - from u=0 to u= 0.125
2nd part: -3/2 ∫ (1 - v^3/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...

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/ '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/ '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/ '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).

 

[ognik]

Hi - I think I have the same bounds, first part with u=t^1/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/2^3/2 which is the same as 2^-3/2. I assume this will make no difference? ... so - now to recheck the program...

 

[ognik]

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
temp=(1.0-(u**3.0))
out=(1/temp)**(1/3.0)
out=3*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
temp=(1.0-(v**3/2.0))
out=(1/temp)**(2/3.0)
out=(3/2)*out
! print 10, 'v= ',v, ' out= ', out !debug
Func_2= out
10 format (2(A,F10.7))
End Function

 

[ognik]

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?

 

[ognik]

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=u³
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-u³)^(-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.

 

[Hawkeye18]

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.

 

[ognik]

Many thanks about the limit, don't 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)

 

[Hawkeye18]

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.

 

[ognik]

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