Evaluating Principal Value Integral

Dear All,

I have a problem to do numerical integration of improper integral. I need to calculate the principal value of integral below, from zero to infinity. My final aim is to evaluate this integral for any values of alpha, H, c, and z. So, the singular point(s) are always changing and it means I can’t just do trial and error to find the limit of the integration to avoid the pole(s), right? The function is depicted below

fun =@(S) (2*(S+alpha).*exp(-S.*H).*cosh(S.*(H+c)).*cosh(S.*(H+z)))./((S.*sinh(S.*H))-(alpha*cosh(S.*H)))
integral(fun,0,inf)

Note: alpha, H, are real positive whereas c, and z are real negative with abs(c) and abs(z) < H

Any idea will be highly appreciated.

Thank you very much before. Regards,

Fredo Ferdian

[fredo ferdian]

Update:

Somehow, I think I can determine the pole location (even I am not really sure whether it works for any arbitrary input H, c and z). For example, I put some value: H = 100 ; alpha = 0.1019 , z = c = -1; and I can determine the pole location by:

 f    = @(S) S.*sinh(S.*H)- alpha.*cosh(S.*H);
 Pole = fzero(f,alpha);

For the integral, I calculate them by:

 fun =@(S)(2.*(S+alpha).*exp(-S.*H).*cosh(S.*(H+c)).*cosh(S.*(H+z)))./((S.*sinh(S.*H))-(alpha.*cosh(S.*H)))
LeftPart   = integral(fun,0,Pole);
RightPart  = integral(fun,(Pole+(realmin)),7);
Results    = LeftPart+RightPart;

The problem here is about the determination of the integral for the RightPart.

1. First problem, I use (Pole+realmin) to exclude the singularity which seems to be too tight tolerance. Maybe I can just increase the tolerance ☹.
2. Second problem, After S = 7 (somewhere around 7.2), my denominator becomes infinite due to the hyperbolic function (because reaching the maximum floating-point value) and subsequently my nominator as well (even not in the same time). This fact gives NaN results since inf/inf = NaN . Therefore, if I use upper limit of integration for the RightPart as inf (instead of 7), the results of the integration will gives NaN. Even actually the integrand is convergence. That is the reason why I put 7 as my upper limit of integration of the RightPart.

Any idea to avoid the inf/inf due to the hyperbolic function? I feel like its kind of pity to loss the accuracy just because of this since the integrand is converge. And considering the fact that I will have various input for H c and z.

Any idea everyone? Thank you very much!

Regards,

Fredo Ferdian

[Torsten]

I did not consider your problem in greater depth, but what makes you think that the principal value of your integral exists and is finite ?

And does the denominator have only one zero ?

Best wishes

Torsten.

[fredo ferdian]

Hi Torsten,

Actually, it is a well defined problem that already solved before. It is clearly stated that I need to take the principal value of that integral. I just want to try my self. About the denominator, actually there is only one real root and infinitely many imaginary roots. The denominator is a dispersion equation of wave if you might already hear before. So I guess, I only need to consider the one real root for the pole right(?). Initially, I just make some upperbound limit where the (S) can capture the maximum floating point before it goes to infinity/NaN. But, apparently after I tried with a lot of different inputs, you are right, after some evaluation, now I hesitate whether the integral is really finite. But I guess this question does not really related with MATLAB it self. Thanks for your sharp reminder anyway.

Answer 1

[David Goodmanson]

Hello fredo,

It makes sense to rescale the independent variable using

s = y/H    ds = dy/H

which leads to

integrand = ((2/H)*(y+aH).*exp(-y).*cosh(y*(1+cH)).*cosh(y*(1+zH))./((y.*sinh(y))-(aH*cosh(y))
where
aH = alpha*H    cH = = c/H    zH = z/H

This effectively reduces the problem from four variables to three since H only enters in as a constant in front. For y somewhat larger than 700 the expression gives NaNs as you mentioned, because the exponentials overflow before division of numerator by denominator. However, multiplying both numerator and denominator by exp(-y) leads to

integrand = (1/H)*(y+aH).*exp(y*(cH+zH)).*(1+exp(-2*y*(1+cH))).*(1+exp(-2*y*(1+zH))) ...
            ./(y.*(1-exp(-2*y))-aH*(1+exp(-2*y)));

This has only decreasing exponentials and takes care of the 'second problem'.

To integrate, here are a couple of methods. Method 1 finds the zero in the denominator and integrates up to a distance y1 on each side. You can vary y1 to see how close you have to get to y0 for the answer to settle down. Method 2 goes around y0 (the pole) with a path in the complex plane, so you don’t need to find the y0 location. The path goes vertically from the origin to (in this case) .1i, and then horizontally from that point out to (inf + .1i). I don’t know for a fact that this path always works but it seems to agree in the cases I tried. (The principal value is the average of two paths: this path which goes above y0 to +inf, and a second path that goes below y0 to +inf. But for a real function these paths are complex conjugates of each other so you end up with the real part of the first path).

In both case the (1/H) factor was ignored in the function and inserted at the end.

H = 100;  alpha = .1019;  c = -1;  z = -1;

% method 1
fun0 = @(y,aH) y.*(1-exp(-2*y))-aH*(1+exp(-2*y));
y0 = fzero(@(x) fun0(x,aH),aH)

y1 = 1e-6;
Ia = integral(@(y) fun(y,aH,cH,zH),0,y0-y1);
Ib = integral(@(y) fun(y,aH,cH,zH),y0+y1,inf);
I1 = (1/H)*(Ia+Ib)
I1 =  0.6296

% method 2
const = .1i;
Ic = integral(@(y) fun(y,aH,cH,zH),0,const);
Id = integral(@(y) fun(y+const,aH,cH,zH),0,inf);
I2 = (1/H)*real(Ic+Id)
I2 = 0.6296

 I1-I2
 ans = -1.0878e-08

function A = fun(y,aH,cH,zH)
  A = (y+aH).*exp(y*(cH+zH)).*(1+exp(-2*y*(1+cH))).*(1+exp(-2*y*(1+zH))) ...
  ./(y.*(1-exp(-2*y))-aH*(1+exp(-2*y)));
end

[fredo ferdian]

Dear David Goodmanson,

Thank you very much for you reply and elegant solution. The code works very well. However, if you don’t mind, I would like to expand the complexity of the problem and understand the behaviour of the function of that integrand. Actually, I have some indicator which can evaluate whether this integrand is calculated correctly or not. Let say, I have problem A which can be solved by 2 solutions:

1. Solution1 = Term1 + Term2 + Term3

2. Solution2 = DifferentFormula

Solution1 and Solution2 supposed to give exactly the same solutions. Solution2 is alternative form of Solution1 which derived by some method so called “Eigenfunction expansion” of Solution1 which I don’t really familiar myself with. PV integral that I asked above, multiplied with TermR is equal to Term3. I successfully created the function of Term1, Term2 and Solution2. TermR is besselj(0,S*R) for different values of R with R>0 and real. So, I slightly modified the integrand which I asked before and implement your solution become:

        aH = alpha*H;
        cH = c/H;
        zH = z/H;
        fun0 = @(y,aH) y.*(1-exp(-2*y))-aH*(1+exp(-2*y));
        y0 = fzero(@(x) fun0(x,aH),aH);
        y1 = 1e-12;

        fun = @(y,aH,cH,zH) (y+aH).*exp(y*(cH+zH)).*(1+exp(-2*y*(1+cH))).*(1+exp(-2*y*(1+zH))) ...
              .*besselj(0,(y./H).*R)./(y.*(1-exp(-2*y))-aH*(1+exp(-2*y)));
        Ia = integral(@(y) fun(y,aH,cH,zH),0,y0-y1);
        Ib = integral(@(y) fun(y,aH,cH,zH),y0+y1,inf);
        I1 = (1/H)*(Ia+Ib);
        Term3     = I1;

I found that in some particular results of z and c, Solution1 and Solution2 give different results at small value of R. Moreover, when I define z=c=0, the solution of Term3 (the PV integrand) will not converge due to the fact that it goes asymptotic to one instead of zero (whereas the upperbound limit of the integration is infinity). So, I tend to conclude that this PV integral evaluation is not really give me correct solution for any arbitrary input of z,c,H, and R (due to the comparison with Solution2). I have two question/asking for some recommendation to approach this matter:

1. I would like to define, where this evaluation is giving me the correct solution. At which particular input of z,c,H and R I can ensure this results.

2. For my final evaluation, I need to do surface integration of above solution. Let say the R is distance of source point and field point. The field point actually is not a point but rather a rectangular panel which can be defined by (u,v). And the solution is the surface integration along this rectangular panel (which certainly has different value of R). So, I wonder, is it possible to solve this by doing something like numerical surface integration of some anonymous function like:

     FinalFunction =@(u,v) Term3(u,v)
     FinalOutput = integral2(FinalFunction,0,1,0,1)

For rectangular panel A bounded from u = 0 to 1 and v bounded from v=0 to 1. Without this surface integration, I can only approximate the final solution by multiply the Solution1 with area of rectangular panel A. I face some difficulties because when defining the Term3(u,v), I need to do integral(@(y) Term3(u,v,y) ) . Is it really possible to do this? I tried to look for several similar questions but most of the solutions suggest to use Symbolic Math Toolbox whereas in this problem, I do not think it is really possible due to the complexity of the PV integral.

Any idea? Or should I ask this by creating new topics?

Regards,

Fredo Ferdian

[fredo ferdian]

To be honest, I dont understand your "second method": the integration with a path in the complex plane. Since the result of integration quite sensitive (with respect to Term1 and Term2), I'm afraid the first method is not really accurate. Even for 10^-3 difference between Method1 and Method2 can gives significant difference between Solution1 and Solution2 (ex, Solution1 with method2 is in a good agreement with respect to Solution2 rather than Solution1 with method1). What is the 0+0.1i represent to? Is there any modification of your code for method2 needed to implement arbitrary input of z,c,H,and R?

Regards,

Fredo Ferdian

[David Goodmanson]

Hi fredo,

In line with what I had before I will set RH = R/H. You went to a new problem when you put the J0 factor into the integrand. However, since the integral from 0 to infinity of J0(y*RH) is bounded, I think that the overall integral converges unless cH=zH=RH=0.

I can't comment on method 1 since none of the details of Solution 1 and Solution 2 are provided. You did not mention why you don't just use Solution 2. (Too slow)?

For method 2, the path of integration starts at the origin and goes up to the point 0+.1i in the complex plane. .1i was chosen somewhat arbitrarily. After that the path of integration goes straight out to the right, to inf+.1i. The second path goes from the origin down to -.1i and then to the right to inf-.1. Complex variable theory says that the answer is the average of these two paths which in this case is the same as the real part of the first path. You will need J0 with complex argument which Matlab does support.

The panel aspect could be difficult unless you can gain a lot of confidence in some method of integration and even then you might have to set up a table of solutions and interpolate or some such.

Here is an example where the answer is known exactly which you could try techniques on:

y(x) = 1/(cosh(x)-cosh(x0))
[obvious 0 in denominator at x = x0]
pv Integral{0,infinity} y(x) dx = - x0/sinh(x0)

[fredo ferdian]

Dear David Goodmanson,

Thank you very much for your suggestion. Yes, I think I will use the method2. I dont think it is possible to provide the code details for Solution1 and Solution2 here, since the formula are quite complicated. Yeah, in the end I use the method 2 as your suggestion. However, I slighly modified the path. I have one last question maybe. Since I know exactly my pole location (k), I made the path just around the pole. For example like this

 % 
tol    = 10^-8;
const1 = k-tol;                                                                
const2 = k+tol*i;                                                                
const3 = k+tol;
const4 = k-tol*i;    
Ic     = integral(@(S) fun(S),0,const1);
Id     = integral(@(S) fun(S),const1,const2);
Ie     = integral(@(S) fun(S),const2,const3);
If     = integral(@(S) fun(S),const4,const3);
Ig     = integral(@(S) fun(S),const1,const4);
Ih     = integral(@(S) fun(S),const3,inf);
I2     = real(Ic+Id+Ie+If+Ig+Ih);

However, the direction of the path does effects the results. I'm sorry for my lackness understanding in evaluating PV integral, but for If and Ig, if I put the integration limit the other way around (eg: If from const3 and const4 and Ig from const4 to const1),the one which compensate each other is the real part instead. Is this correct? I feel weird because I thought I only need to define close boundary around the pole.

Regards,

Fredo Ferdian

[David Goodmanson]

Hello Fredo,

You have a reasonable approach except that you don’t have to have tol anywhere near as small as 1e-8. In fact that might be counterproductive.

Before, I used y as the variable of integration. Pictorially it makes more sense to use x, so that's done here.

For the p.v. you integrate on one path that goes over the pole and on another path that goes under, and take the average. You have been looking at 'rectangular' paths that would look like

                    d------e  
  path1    a=0------b  x0  f-------g=inf      pole is at x0 (called k before)
.
  path2    a=0------b  x0  f-------g=inf
                    m------n

Paths ab and fg on the x axis give real integrals, and the 'curved' paths bdef and bmnf give complex integrals.

As I mentioned you don't have to actually average the two paths, since taking the real part of either path by itself gives the same answer.

For the points on either side of x0, let b = x0-r, f = x0+r, where notation r is used instead of tol. If r is really small like 1e-8, then the integrand takes on huge values near the pole, on the order of 1e8. That's not a real bad problem with an adaptive integration routine like 'integral', but the highly varying integrand affects the accuracy and also takes more time. And there is no reason to do that, because after including the curve integrals the answer is just as exact when, for example r = .1 .

Instead of a three-part rectangular path it's easier to integrate along a semicircular path of radius r, centered at x0, going from b to f. It's an integral in variable theta, with

z = x0 + r*exp(i*theta)        
dz = (dz/dtheta)*dtheta   where
dz/dtheta = i*r*exp(i*theta)             note the important factor of i.

Theta is measured counterclockwise from the x axis as usual and goes from pi to 0 for the upper curve and -pi to 0 for the lower curve.

The code below is similar to before, except the bessel function factor is included. Method 3 is the path discussed above, and method 2 is as before with the path

          d-----------------------d+inf  
         a=0         x0

You can independently adjust r in method 3, and d in method 2, to see how the answer is affected. Over a wide range from 1e-1 to 1e-8 or less for each of them, the two methods agree pretty well and the computations for the larger values of r are faster.

H = 100;  alpha = .1019;  c = -1;  z = -1;  R = 10;
aH = alpha*H;    cH = c/H;    zH = z/H;    RH = R/H;

% find pole
fun0 = @(x,aH) x.*(1-exp(-2*x))-aH*(1+exp(-2*x));  % the denominator
x0 = fzero(@(x) fun0(x,aH),aH);

% method 2 (from before) path in complex plane does not need pole location
d = 1e-1;
% integrate straight up from 0, then across to the right
Iup = integral(@(x) funJ(x,aH,cH,zH,RH),0,i*d);
Iright = integral(@(x) funJ(x+i*d,aH,cH,zH,RH),0,inf);
I2 = (1/H)*real(Iup+Iright);

% method 3 (integrate to each side of pole, then around the top in a semicircle
r = 1e-1;
Iab = integral(@(x) funJ(x,aH,cH,zH,RH), 0, x0-r);
Ieh = integral(@(x) funJ(x,aH,cH,zH,RH), x0+r, inf);
Iabove = integral(@(theta) funJarc(theta,r,x0,aH,cH,zH,RH), pi, 0);
I3 = (1/H)*(Iab + Ieh + real(Iabove));

I2and3 = [I2 I3]
diff23 = I2-I3

function A = funJ(x,aH,cH,zH,RH)
    A = (x+aH).*besselj(0,x*RH).*exp(x*(cH+zH)) ...
      .*(1+exp(-2*x*(1+cH))).*(1+exp(-2*x*(1+zH))) ...
      ./(x.*(1-exp(-2*x))-aH*(1+exp(-2*x)));
end

function A = funJarc(theta,r,x0,aH,cH,zH,RH)
    zcirc = x0 + r*exp(i*theta);
    % dz/dtheta = i*r*exp(i*theta);
    A = (i*r*exp(i*theta)).*funJ(zcirc,aH,cH,zH,RH);
end