help-octave
[Top][All Lists]
Advanced

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

Re: Double integration results in "error: quad: invalid recursive call"


From: Geordie McBain
Subject: Re: Double integration results in "error: quad: invalid recursive call"
Date: Fri, 25 Feb 2005 09:44:08 -0500
User-agent: Mutt/1.5.6+20040907i

On Thu, Feb 24, 2005 at 01:46:11AM -0500, address@hidden wrote:
> according the method you use in colloc() in last email, I change my variable 
> x, y to rsin(theta), rcos(theta), it indeed make problem become 2 independent 
> variable if I chose to integrate on a circle (area not perimeter) radius 3, 
> so 
> f= 9x^2 -3z
> 
> ---------------------------------------------------
> octave:1> function z=f(x,y)
> >                  z= (9* x^2 * (sin(y)^2) - 3 * x * cos(y)) * x;
> > endfunction
> octave:2> n=1; [r,A,B,q]=colloc(n); r*=3; r+=0; q*=2*pi; q'*f(r,r')*q
> ans = 1174.3
> octave:3> 
> -----------------------------------------------------------
> but that is not the answer, which should be 572.555
> 
> any comment?
> eric

Yes, the method I demonstrated was concise but doesn't generalize
quite so simply.  It goes like this.  The basic thing you need to know
is how the colloc n, r, and q work for approximating integrals:

  integral from x=0 to x=1 of f(x) ~= sum over i=1:n of q(i)*f(r(i))

In Octave, the sum on the RHS can be written as q'*f(r), provided your
function f returns a column vector the same size as r.

For an arbitrary interval, you need to (i) map the nodes and (ii)
correspondingly multiply the weights:

  integral_{x=a}^{x=b} f(x) ~= sum_{j=1}^n {(b-a)*q(j)} * {f((b-a)*r(j)+a)}

to do an integral over a rectangle, you need to stretch in both directions, so
 
  int_{x=a}^b int_{y=c}^d f(x,y) ~= 

        sum_{i=1}^n sum_{j=1}^n {

                {(b-a)*q(j)} * {{c-d)*q(i)} * f((b-a)*r(j)+a, (c-d)*r(i)+c)

        }

Provided the function f(x, y) returns a matrix with
z(i,j)=f(x(i),y(j)), you can write the double sum as

        (b-a)*q' * f((b-a)*r+a, (c-d)*r+c) * (c-d)*q

However, your function f doesn't do that: check it's individual entries.

Basically, the things you need to fix are:

        1) map each independent variable independently

        2) only use q'*f*q form if f is the correct matrix

Also, as I noted in the original message, you need to increase n, the
number of nodes in each direction.

I run your problem like this: define f:

function z = f (r, theta)
  z = 9 * r.^3 * (cos (theta)).^2 - 3 * r.^2 * sin (theta);
endfunction

then:

  octave> n=1; [t,A,B,q]=colloc(n); 3*q'*f(3*t,2*pi*t')*2*pi*q
  ans = 572.56

Looks good, right?  But don't be fooled, since:

  octave> n=2; [t,A,B,q]=colloc(n); 3*q'*f(3*t,2*pi*t')*2*pi*q
  ans = 66.299

Don't worry: the method does converge:

 octave> for n=1:20; [t,A,B,q]=colloc(n); disp ([n, 
3*q'*f(3*t,2*pi*t')*2*pi*q]); endfor
    1.0000  572.5553
   2.0000  66.2988
    3.0000  875.9876
    4.0000  500.5666
    5.0000  582.0991
    6.0000  571.7409
    7.0000  572.6039
    8.0000  572.5531
    9.0000  572.5553
   10.000  572.555
   11.000  572.555
   12.000  572.555
   13.000  572.555
   14.000  572.555
   15.000  572.555
   16.000  572.555
   17.000  572.555
   18.000  572.555
   19.000  572.555
   20.000  572.555

O.K.?  

This is quite a low-level way to do integration in Octave, so maybe
it's not what you're used to or what you're after; however, it is
powerful and efficient if you're careful.

Geordie McBain
www.aeromech.usyd.edu.au/~mcbain



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------



reply via email to

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