Re: complex line integral in octave

[Urs Hackstein]

Thanks a lot for your explanations, they were very helpful. I have still a problem concerning the numerical integration. I want to do the computation for several different polynomials f, so I have seven parameters La,Lb,Lc,Ia,Ib,Ic, RC besides the integration variable t. My attempt for the integral over the first edge of the rectangle was the following:

function g1=complexfunction1(La,Lb,Lc,Ia,Ib,Ic,RC,t)
sd=z(1,-10.^11,1,poly(La,Lb,Lc,Ia,Ib,Ic,RC));
wmax=3*10.^11;
a0=koeff0(La,Lb,Lc,Ia,Ib,Ic,RC);
a1=koeff1(La,Lb,Lc,Ia,Ib,Ic,RC);
a2=koeff2(La,Lb,Lc,Ia,Ib,Ic,RC);
a3=koeff3(La,Lb,Lc,Ia,Ib,Ic,RC);
a4=koeff4(La,Lb,Lc,Ia,Ib,Ic,RC);
a5=koeff5(La,Lb,Lc,Ia,Ib,Ic,RC);
a6=koeff6(La,Lb,Lc,Ia,Ib,Ic,RC);
G1=ones(1,6).*(sd-10.^8+2*10.^8*t).^(1:6);
H1=[a1 2*a2 3*a3 4*a4 5*a5 6*a6];
G11=ones(1,7).*(sd-10.^8+2*10.^8*t).^(0:6);
H11=[a0 a1 a2 a3 a4 a5 a6]
g1=(7*(sd-10.^8+2*10.^8*t).^7+H1.*G1')*2*10.^8/((sd-10.^8+2*10.^8*t).^7+H11*G11')
endfunction

 Here sd is the real part of the root inside the rectangle and the ai are the coefficients of the polynomial.

function g1r=real_part1(La,Lb,Lc,Ia,Ib,Ic,RC,t)
g1r=real(complexfunction1(La,Lb,Lc,Ia,Ib,IC,RC,t))
endfunction

function g1i=imag_part1(La,Lb,Lc,Ia,Ib,Ic,RC,t)
g1i=imag(complexfunction1(La,Lb,Lc,Ia,Ib,IC,RC,t))
endfunction

function complex_answer1=complex1(La,Lb,Lc,Ia,Ib,Ic,RC)
real_integral1=quad("real_part1",0,1)
imag_integral1=quad("imag_part1",0,1)
complex_answer1=real_integral1+I*imag_integral1
endfunction

I guess that quad doesn't know which variable is the integration variable, but how can we fix it?

Thanks a lot in advance.

2013/4/12 Stephen Montgomery-Smith <address@hidden>

On 04/12/2013 05:29 AM, Urs Hackstein wrote:
> Dear Stephen,
>
> thanks a lot for your long reply. Unfortunately, we are talking about
> different integrals. Your are dealing with the zero-counting integral
> \int_\alpha f ' (s)/f(s) ds, I work on the integral \int_\alpha s*f '
> (s)/f(s) ds, where * is the common multiplication. From complex function
> theory we know that it equals a, if a is the only root of f in the
> rectangle. We want to compute this root, thus method 2 doesn't work for
> us. Remains method 1: what is the best (most efficient and stable) built
> in integration command?
>

First, you meant 2 pi I a, not a.  I'm sure it was a typo, but it is the
kind of typo that propagates into code and creates a hard to trace bug.
 (At least, that is what it does for me.)

Second, I think you have a fairly straightforward integral to compute,
so I think any numerical method should do fine.  (Unless one of the
roots lies very close or is on the curve.)

Third, if you know that there is only one root inside the rectangle, you
could use the output of Method 1 as the input to Newton's Method, which
will very quickly and accurately converge to the desired root.

Finally, friends of mine highly recommend the NIntegrate command in
recent versions of Mathematica (version 8 or higher).  It might even
handle the case where one root lies on the curve (where presumably it
correctly computes the principle value).  So perhaps you could check
your answers using the wolfram alpha web site,

2013/4/12 Stephen Montgomery-Smith <address@hidden>

On 04/12/2013 05:29 AM, Urs Hackstein wrote:
> Dear Stephen,
>
> thanks a lot for your long reply. Unfortunately, we are talking about
> different integrals. Your are dealing with the zero-counting integral
> \int_\alpha f ' (s)/f(s) ds, I work on the integral \int_\alpha s*f '
> (s)/f(s) ds, where * is the common multiplication. From complex function
> theory we know that it equals a, if a is the only root of f in the
> rectangle. We want to compute this root, thus method 2 doesn't work for
> us. Remains method 1: what is the best (most efficient and stable) built
> in integration command?
>

First, you meant 2 pi I a, not a.  I'm sure it was a typo, but it is the
kind of typo that propagates into code and creates a hard to trace bug.
 (At least, that is what it does for me.)

Second, I think you have a fairly straightforward integral to compute,
so I think any numerical method should do fine.  (Unless one of the
roots lies very close or is on the curve.)

Third, if you know that there is only one root inside the rectangle, you
could use the output of Method 1 as the input to Newton's Method, which
will very quickly and accurately converge to the desired root.

Finally, friends of mine highly recommend the NIntegrate command in
recent versions of Mathematica (version 8 or higher).  It might even
handle the case where one root lies on the curve (where presumably it
correctly computes the principle value).  So perhaps you could check
your answers using the wolfram alpha web site,

2013/4/12 Stephen Montgomery-Smith <address@hidden>

On 04/12/2013 05:29 AM, Urs Hackstein wrote:
> Dear Stephen,
>
> thanks a lot for your long reply. Unfortunately, we are talking about
> different integrals. Your are dealing with the zero-counting integral
> \int_\alpha f ' (s)/f(s) ds, I work on the integral \int_\alpha s*f '
> (s)/f(s) ds, where * is the common multiplication. From complex function
> theory we know that it equals a, if a is the only root of f in the
> rectangle. We want to compute this root, thus method 2 doesn't work for
> us. Remains method 1: what is the best (most efficient and stable) built
> in integration command?
>

First, you meant 2 pi I a, not a.  I'm sure it was a typo, but it is the
kind of typo that propagates into code and creates a hard to trace bug.
 (At least, that is what it does for me.)

Second, I think you have a fairly straightforward integral to compute,
so I think any numerical method should do fine.  (Unless one of the
roots lies very close or is on the curve.)

Third, if you know that there is only one root inside the rectangle, you
could use the output of Method 1 as the input to Newton's Method, which
will very quickly and accurately converge to the desired root.

Finally, friends of mine highly recommend the NIntegrate command in
recent versions of Mathematica (version 8 or higher).  It might even
handle the case where one root lies on the curve (where presumably it
correctly computes the principle value).  So perhaps you could check
your answers using the wolfram alpha web site,