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*).