In the example examples/quad/TjTkDisk.html
, we illustrated formulas due to R. M. Slevinsky for the integral over the unit disk of a project of Chebyshev polynomials $T_j(x) T_k(y)$. The fascinating property of such integrals is that they are always equal to zero except when $j$ and $k$ and both even and differ by $0$ or $\pm 2$.
We also commented at the end of that example that these formulas could be used as the basis of a Chebfun2 command sumdisk
, which would elegantly compute the double integral of a chebfun2, not over its square (or rectangle) of definition but over the inscribed disk (or elliptical region). Subsequently, such a code has been written by the first author. Here we show it off.
For a trivial example, suppose our chebfun2 is the constant 1. Its integral over the square is 4,
f = chebfun2(@(x,y) 1); format long sum2(f)
ans = 4
but its integral over the disk is $\pi$,
sumdisk(f)
ans = 3.141592653589793
As another example, let's consider the bivariate Gaussian $\exp(-(x^2+y^2)/2$. Here is its integral over the unit disk:
cheb.xy f = exp(-(x.^2+y.^2)/2); sumdisk(f)
ans = 2.472240777719226
Switching to polar coordinates enables us to perform the integral exactly; it is $2\pi(1 - 1/\sqrt e)$:
exact = 2*pi*(1-exp(-.5))
exact = 2.472240777719227
We must make a comment about the significance of sumdisk
. We would certainly not claim that a competitive way to integrate a function over a disk is to make a chebfun2 of it and then call sumdisk
. It would be much more efficient to work on the integral over the disk directly, and indeed, examples/quad/TjTkDisk.html
gives a sample code for doing just that, which we illustrate again here:
fpolar = @(r,t) f(r*cos(t),r*sin(t)); fr = @(r) r*sum(chebfun(@(t) fpolar(r,t),[0,2*pi],'trig')); I = sum(chebfun(@(r) fr(r),[0 1]))
I = 2.472240777719227
The point of sumdisk
is two-fold: it shows off some elegant mathematics, and it provides an good way to compute integrals over a disk if, for whatever reason, you are already working with chebfun2 objects on a square.
Here's another example. Suppose $f$ is a harmonic function, which for convenience we might obtain as the real part of an analytic function. Chebfun2 can do this very conveniently, like this
fcomplex = chebfun2(@(z) cos(2*cosh(z))); f = real(fcomplex); plot(f), camlight
Here we use sumdisk to compute the mean of $f$ over the unit disk:
sumdisk(f)/pi
ans = -0.416146836547142
Since $f$ is harmonic, this must be the same as the value of $f$ at the origin:
f(0,0)
ans = -0.416146836547142
Another option for such integrals is to use Diskfun (see chapter 16 of the Chebfun Guide).