### 1. Numerical integration over the disk

In studying cubature formulas, we needed to compute the integrals of products of Chebyshev polynomials $T_j(x) T_k(y)$ over the unit disk, like this one:

T8 = chebpoly(8); T16 = chebpoly(16);
s = linspace(-1,1,160); [xx,yy] = meshgrid(s,s);
ff = T8(xx).*T16(yy);
ff(xx.^2+yy.^2>1) = NaN;
contour(s,s,ff), axis equal off
hold on, plot(chebfun('exp(1i*x)',[0 2*pi]),'k'), hold off If $j$ or $k$ is odd, the integral is zero, but even if they are both even, to our surprise, we found that the integrals are still usually zero! In fact, a nonzero result only shows up if $j$ and $k$ differ by 0 or 2. Thus the function plotted above, for example, has integral exactly zero over the disk. This is not obvious.

Speaking in general, suppose we want to integrate a smooth function $f(r,t)$ numerically over the unit disk, where $r$ is radius and $t$ is angle. One could use Diskfun, but here, we use standard Chebfun. Let's assume that $f$ is defined by a function that accepts a scalar $r$ and a vector $t$, like this one:

f = @(r,t) 1 + 0*t;

Here is a numerical confirmation that the integral of $1$ is $\pi$:

format long
fr = @(r) r*sum(chebfun(@(t) f(r,t),[0,2*pi],'trig'));
I = sum(chebfun(@(r) fr(r),[0 1]))
Iexact = pi
I =
3.141592653589793
Iexact =
3.141592653589793


For the function $f(r,t) = r^2$, the integral is $\pi/2$:

f = @(r,t) r^2 + 0*t;
fr = @(r) r*sum(chebfun(@(t) f(r,t),[0,2*pi],'trig'));
I = sum(chebfun(@(r) fr(r),[0 1]))
Iexact = pi/2
I =
1.570796326794897
Iexact =
1.570796326794897


For the function $f(r,t) = r^2 \cos^2(t)$, the integral is $\pi/2$:

f = @(r,t) r^2*cos(t).^2;
fr = @(r) r*sum(chebfun(@(t) f(r,t),[0,2*pi],'trig'));
I = sum(chebfun(@(r) fr(r),[0 1]))
Iexact = pi/4
I =
0.785398163397448
Iexact =
0.785398163397448


### 2. Numerical integration of products of Chebyshev polynomials

What about those products of Chebyshev polynomials? Here is a matrix showing the numerically computed integrals for $k = 0,2,4,6,8,10$. As claimed above, the matrix is tridiagonal.

tic
I = zeros(6);
format short
for j = 0:2:10
Tj = chebpoly(j);
for k = 0:2:j
Tk = chebpoly(k);
f = @(r,t) Tk(r*cos(t)).*Tj(r*sin(t));
fr = @(r) r*sum(chebfun(@(t) f(r,t),[0,2*pi],'trig'));
I(1+j/2,1+k/2) = sum(chebfun(@(r) fr(r),[0 1]));
end
end
I = I + tril(I,-1)'
time_elapsed_in_seconds = toc
I =
3.1416   -1.5708   -0.0000    0.0000   -0.0000    0.0000
-1.5708    0.5236    0.2618    0.0000   -0.0000   -0.0000
-0.0000    0.2618   -0.1047   -0.1571    0.0000   -0.0000
0.0000    0.0000   -0.1571    0.0449    0.1122    0.0000
-0.0000   -0.0000    0.0000    0.1122   -0.0249   -0.0873
0.0000   -0.0000   -0.0000    0.0000   -0.0873    0.0159
time_elapsed_in_seconds =
3.7174


### 3. Analytic expressions for the integrals

Let $I_{jk}$ denote the integral of $T_j(x) T_k(y)$ over the unit disk. The following explicit formulas derived by the first author (details not reported here) give the integrals:

$$I_{00} = \pi,$$

$$I_{02} = I_{20} = -{\pi\over 2},$$

$$I_{kk} = {\pi (-1)^{k/2} \over 2 - 2k^2} \quad \hbox{(k even, k\ge 2)},$$

$$I_{k,k+2} = I_{k+2,k} = {\pi (-1)^{1+k/2} \over 4k + 4} \quad \hbox{(k even, k\ge 2)},$$

In all other cases $I_{jk} = 0$.

Using these formulas, we can reproduce the matrix above as follows:

I = zeros(6);
I(1,1) = pi;
I(2,1) = -pi/2;
for k = 2:2:10
I(1+k/2,1+k/2) = pi*(-1)^(k/2)/(2-2*k^2);
end
for k = 2:2:8
I(2+k/2,1+k/2) = -pi*(-1)^(k/2)/(4*k+4);
end
I = I + tril(I,-1)'
I =
3.1416   -1.5708         0         0         0         0
-1.5708    0.5236    0.2618         0         0         0
0    0.2618   -0.1047   -0.1571         0         0
0         0   -0.1571    0.0449    0.1122         0
0         0         0    0.1122   -0.0249   -0.0873
0         0         0         0   -0.0873    0.0159


### 4. Application to integration of general functions

The results of the last section imply that there is an immediate way to compute the integral of a chebfun2 over a disk: just take the appropriate linear combination of its bivariate Chebyshev coefficients. After this example was initially written, we developed this idea into a Chebfun2 sumdisk command. See the example "Sumdisk for integration over a disk".