Re: complex line integral in octave

[Stephen Montgomery-Smith]

This is a programming issue, not mathematics.  If you look at "help
quad" you will see that the function takes only one parameter.  The
other parameters, La,Lb,Lc,Ia,Ib,Ic,RC, will have to be passed in some
other manner.  My suggestion would be to use global variables.  But
others more proficient in octave or matlab may have better suggestions.

Also, if you want to use other integration algorithms, some of them
require the function to accept vectors, and then evaluate the function
at all those values.  And you may as well try them out to see if they
give better or faster results.

Also octave can handle complex numbers, so I wouldn't mess with
splitting up the real and imaginary parts.

To compute the line integral from z=a to z=b, I would use the following
code.  But it isn't tested.  I also pass a and b using global variables.

function g = myf(t)
global La,Lb,Lc,Ia,Ib,Ic,RC
global a,b

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

z = a + (b-a)*t;

f = (((((a6.*z.+a5).*z.+a4).+a3).*z.+a2).*z.+a1).*z.+a0;
fp = ((((6*a6.*z.+5*a5).*z.+4*a4).+3*a3).*z.+2*a2).*z.+a1;

g = z.*fp./f

end

On 04/15/2013 08:51 AM, Urs Hackstein wrote:
> 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
> <mailto: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
> <mailto: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
> <mailto: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,
> 
>