[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Using quad() for multidimensional integration
From: |
Ted Harding |
Subject: |
Re: Using quad() for multidimensional integration |
Date: |
Thu, 17 Aug 1995 19:22:25 +0200 (BST) |
( Re Message From: Vinayak Dutt )
>
> I just found that even that simple double integral does not work OK.
>
> But I still wonder why would quad run out of stack space on simple
> double integrals (forget triple integral). I don't see that many levels
> of recursion here?
>
Nor do I. But they happen. Or at least there is a totally unnecessary
number of calls to the function being integrated. Take the simple example
function u=U(x)
fprintf(stdout,"U(%4.2f)\n",x);
u=0.5*x*x;
endfunction
and do " quad('U',0,1) " and you get (packed onto three lines):
U(0.50) U(0.01) U(0.99) U(0.07) U(0.93) U(0.16) U(0.84) U(0.28) U(0.72)
U(0.43) U(0.57) U(0.00) U(1.00) U(0.03) U(0.97) U(0.11) U(0.89) U(0.22)
U(0.78) U(0.35) U(0.65) ans = 0.16667
i.e. 21 calls in order to integrate a quadratic! Now (if octave used, like
matlab, Simpson's rule) Simpson's rule is theoretically EXACT IN ONE PASS
for a quadratic. That is to say that if
f(x) = a + b*x + c*x^2
then \int_0^1 f(x) dx = [f(0) + 4*f(0.5) + f(1)]/6 exactly.
Of course "quad" wouldn't know it was quadratic, but on splitting the
range into 2 halves (which requires 2 further calls to U) and combining
the integrals (also exact) over the two halves, and comparing the second
result with the first, convergence would be exact. So just 5 calls are
required.
[NOTE: when I run this in matlab, U is evaluated at 4 points only, namely
0, 0.25, 0.125 and 0.625. Hmm]
Octave "help -i quad" says the integration is done using "Quadpack".
I have a feeling I would like to know more about this "Quadpack".
I don't have the octave source (I downloaded the binary only).
Does anyone out there know how octave does "quad"?
Ted. (address@hidden)