[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
-------------------------------------------------------------