help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: complex line integral in octave


From: Urs Hackstein
Subject: Re: complex line integral in octave
Date: Tue, 16 Apr 2013 15:35:47 +0200

There was another typo, since I guess you meant
f = (((((a6.*z.+a5).*z.+a4).*z.+a3).*z.+a2).*z.+a1).*z.+a0;
fp = ((((6*a6.*z.+5*a5).*z.+4*a4).*z.+3*a3).*z.+2*a2).*z.+a1;

g = z.*(fp.+7.*z.^6)./(f.+z.^7) .*(b-a)

I have still a programming problem. Since I want to call the value of the integral for variable La, Lb,Lc,Ia,Ib,Ic,RC from another greater program, I need something like
function intsteuer1=integralsteuer1(La,Lb,Lc,Ia,Ib,Ic,RC)
as a function that gives back the value of the integral.

If I do

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

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 = d1 + (d2-d1)*t;

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

g = z.*(fp.+7.*z.^6)./(f.+z.^7) .*(d2-d1)

endfunction

and

function intsteuer1=integralsteuer1(La,Lb,Lc,Ia,Ib,Ic,RC)
global La,Lb,Lc,Ia,Ib,Ic,RC
sd=z(1,-10.^11,1,poly(La,Lb,Lc,Ia,Ib,Ic,RC))
global d1=sd-10.^8
global d2=sd+10.^8
intsteuer1=quad("myf",0,1)
endfunction,

i get the following error:

error: for A^b, A must be square
error: called from:
error:   /home/hackstein/Dominantepole5/koeff6.m at line 2, column 7
error:   /home/hackstein/Dominantepole5/poly.m at line 3, column 5
error: evaluating argument list element number 4
error: evaluating argument list element number 1
error:   /home/hackstein/Dominantepole5/integralsteuer1.m at line 3, column 3

Could you give me a hint what I am doing wrong? (sd is a real number, the real part of the only zero in the rectangle.)

2013/4/16 Stephen Montgomery-Smith <address@hidden>
On 04/15/2013 05:02 PM, Stephen Montgomery-Smith wrote:
> 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

I mean

g = z.*fp./f.*(b-a)

>
> 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,
>>
>>
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave
>
>

_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave



2013/4/16 Stephen Montgomery-Smith <address@hidden>
On 04/15/2013 05:02 PM, Stephen Montgomery-Smith wrote:
> 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

I mean

g = z.*fp./f.*(b-a)

>
> 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,
>>
>>
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave
>
>

_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave



2013/4/16 Stephen Montgomery-Smith <address@hidden>
On 04/15/2013 05:02 PM, Stephen Montgomery-Smith wrote:
> 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

I mean

g = z.*fp./f.*(b-a)

>
> 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,
>>
>>
>
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://mailman.cae.wisc.edu/listinfo/help-octave
>
>

_______________________________________________
Help-octave mailing list
address@hidden
https://mailman.cae.wisc.edu/listinfo/help-octave


reply via email to

[Prev in Thread] Current Thread [Next in Thread]