Next: The Eigenvalue Problem Up: Spectral Methods Previous: Spectral Methods

A Diffusion Problem

• Consider the following slightly more complex diffusion equation:

 (3.24)

The equation is more complicated because the diffusion coefficient D is a function of x and y.
• Assume that we're simulating heat flow in a beam of infinite length and a rectangular cross section with sides lx and ly.
• We can always   change our units of length and time so that the beam cross section becomes a square with a side length of - this will lead to certain FFT conveniences. The transformation of coordinates is:
 x' = (3.25) y' = (3.26)

• This implies
 = (3.27) = (3.28)

• Substitute into equation (3.24) to get:
 = = (3.29)

where lr = lx/ly.
• The Diffusion coefficient   has units of . In the above equations we have already replaced x with x' and y with y', but D is still expressed in terms of x and y. Expressing D in terms of x' and y' will throw out some factor, say, Dr and D itself will become D'. In summary D = Dr D', where Dr is measured in and:

 (3.30)

• The factor can be now absorbed into a new unit of time (observe that of Dr cancels with lx2 so that the new unit of time still has the corect dimension of some fraction of seconds):

 (3.31)

so that

 (3.32)

• Changing the unit of time will scale additionally D'. This time we'll absorb that change into a new D' so that no further modifications of the equation itself are needed.
• This restructured equation now has that nice property that whatever the original dimensions of a rectangular cross section of the beam, in the new units the beam extends from 0 to both in x' and in y'. Furthermore distances in x' and in y' are dimensionless.
From this point onwards we drop the primes!
• Such rewriting of PDEs in computationally more convenient although not always natural units is very common. It makes equations more tractable numerically.
• Assume the following initial condition:

 T(x, y, 0) = T0(x, y) (3.33)

• Assume the following fixed boundary conditions:

 (3.34)

• Given our boundary conditions T(x, y, t) can be always represented by:

 (3.35)

because for x=0 or and for y=0 or too. The factor is there for future convenience. You can always shove it into ak'j' if you don't like it.
• Given T(x, y, 0) = T0(x, y) it is easy to evaluate ak'j'(0). First observe that

 (3.36)

You can check this result easily with Emacs Calc:
mr
'integ(sin(2 x) sin(5 x), x, 0, pi)
1:  0
'integ(sin(4 x) sin(4 x), x, 0, pi)
1:  0.5 pi

You can try various other combinations of n and m, but be warned that Emacs Calc may take a long time to do the integral for high values of n and/or m. You can also ask Emacs Calc to evaluate the integral for some general n and m:
'integ(sin(n x) sin(m x), x, 0, pi)
1:  (n cos(n pi) sin(m pi) / m^2
- cos(m pi) sin(n pi) / m)
/ (1 - n^2 / m^2)
'integ(sin(n x) sin(n x), x, 0, pi)
1:  (n pi - cos(n pi) sin(n pi)) / 2 n

It is now clear that for an integer value of n and m the first integral must vanish, and the second integral always returns .

We make use of the above as follows:

 (3.37)

Now you can also see why we have defined our solution with that factor.

• What we've done so far is not Fourier Transform yet. It's just a decomposition into sines.
• Before we switch to a fully blown Fourier Transform let us interrogate Emacs Calc about a few more integrals.
mr
'integ(cos(a x) sin(b x), x, 0, 2 pi)
1:  (cos(2 a pi) cos(2 b pi) / b
+ a sin(2 a pi) sin(2 b pi) / b^2)
/ (a^2 / b^2 - 1) + 1 / b*(1 - a^2 / b^2)

If a and b are two different integers, m and n, then

 (3.38)

for m = n we get a different formula:
'integ(cos(a x) sin(a x), x, 0, 2 pi)
1:  sin(2 a pi)^2 / 2 a

If a is an integer then this is 0, hence:

 (3.39)

Another integral we need to evaluate is
'integ(sin(a x) sin(b x), x, 0, 2 pi)
1:  (b cos(2 b pi) sin(2 a pi) / a^2
- cos(2 a pi) sin(2 b pi) / a)
/ (1 - b^2 / a^2)

This is clearly 0 for integer values of a and b:

 (3.40)

And finally:
'integ(sin(a x) sin(a x), x, 0, 2 pi)
1:  (2 a pi - cos(2 a pi) sin(2 a pi)) / 2 a

which for an integer value of a is simply :

 (3.41)

• We can now extend T0(x, y) analytically to a region by keeping the same coefficients, but letting x and y run between 0 and each.
• Now consider the following integral:

This will produce a variety of terms. First there will be terms containing or and they will all vanish on integration over , including two imaginary terms. The only ones that will survive the slaughter will be the terms with sine functions only, i.e., the ones proportional to :
 (3.42)

• Let us now plug T(x, y, t) as given by quation (3.37) into our diffusion equation (3.34):

where .
• The next step is to multiply both sides by and integrate over . This procedure really is the same as acting with a basis onto a vector in order to obtain a selected vector coefficient. Let us do the left hand side first:

The right hand side is a little more complicated because of the presence of differential operators and an unknown function D(x,y):

• Before going any further let us evaluate

first. It is easy to see that this is:

• Now plug this into the integral to get this very lengthy expression:

• We can't differentiate D because we don't know what it is, but we can integrate by parts and transfer and onto . The result is:
 (3.43)

• Recall that when we had a similar formula for T0(x, y) we have switched to a Fourier transform over . As it turned out then the only difference was that the integral over was twice the integral over . This is also going to be the case here, if D(x,y) is defined to be even about and we can always do that. So in preparation to a move to a Fourier transform we'll extend the integration to the whole . In summary the factor in front of the whole expression will be replaced with , because the 4 will become absorbed into the extended integrals:
 (3.44)

• The next step is to convert equation (3.46), which is still expressed in terms of sines and cosines into a Fourier transform expression.

• The conversion is quite tedious and we will help ourselves with Calc in this process.
• First type into Calc the formula enclosed by large round brackets, i.e., all that's under the integral:
'k k' D cos(k x) cos(k' x) sin(j y) sin(j' y) \
+ l j j' D cos(j y) cos(j' y) sin(k x) sin(k' x)
1:  k k' D cos(k x) cos(k' x) sin(j y) sin(j' y)
+ l j j' D cos(j y) cos(j' y) sin(k x) sin(k' x)

• Next memorize it with ss:
ss
Store: var-F

• Now we will apply the following transformations to this expression:

ar
Rewrite rule(s): cos(a x) cos(b x) := ( cos((a - b) x) + cos((a + b) x) ) / 2,\
sin(a x) sin(b x) := ( cos((a - b) x) - cos((a + b) x) ) / 2
Working...
1:  (cos((j - j') y) - cos((j + j') y)) k
*(cos((k - k') x) + cos((k + k') x)) D k' / 4
+ (cos((k - k') x) - cos((k + k') x)) j
*(cos((j - j') y) + cos((j + j') y)) D j' l / 4

• The next step is to convert these cosines into exponentials using:

ar
Rewrite rule(s): cos(a x) := (exp(i a x) + exp(-i a x)) / 2
Working...
1:  ((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
- (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2) k
*((exp(i x*(k - k')) + exp(-(i x*(k - k')))) / 2
+ (exp(i x*(k + k')) + exp(-(i x*(k + k')))) / 2) D k' / 4
+ ((exp(i x*(k - k')) + exp(-(i x*(k - k')))) / 2
- (exp(i x*(k + k')) + exp(-(i x*(k + k')))) / 2) j
*((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
+ (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2) D j' l / 4
ss
Store: var-F2
as
Working...
1:  D k k'
*((exp(i x*(k - k')) + exp(-(i x*(k - k')))) / 2
+ (exp(i x*(k + k')) + exp(-(i x*(k + k')))) / 2)
((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
- (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2) / 4
+ D j j' l
*((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
+ (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2)
((exp(i x*(k - k')) + exp(-(i x*(k - k')))) / 2
- (exp(i x*(k + k')) + exp(-(i x*(k + k')))) / 2) / 4
ss
Store: var-F3

• Now we'll introduce the following substitutions

The Calc command for doing this type of substitutions is ab or subst. The substitutions are purely structural, i.e., the substituted objects must match verbatim, so this is not like a rewriting rule. Here's how it works on our example:
ab
Substitute old: k - k', new: k_m
Working...
1:  D k k'
*((exp(i x k_m) + exp(-(i x k_m))) / 2
+ (exp(i x*(k + k')) + exp(-(i x*(k + k')))) / 2)
((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
- (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2) / 4
+ D j j' l
*((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
+ (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2)
((exp(i x k_m) + exp(-(i x k_m))) / 2
- (exp(i x*(k + k')) + exp(-(i x*(k + k')))) / 2)
/ 4
ab
Substitute old: k + k', new: k_p
Working...
1:  D k k'
*((exp(i x k_m) + exp(-(i x k_m))) / 2
+ (exp(i x k_p) + exp(-(i x k_p))) / 2)
((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
- (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2) / 4
+ D j j' l
*((exp(i y*(j - j')) + exp(-(i y*(j - j')))) / 2
+ (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2)
((exp(i x k_m) + exp(-(i x k_m))) / 2
- (exp(i x k_p) + exp(-(i x k_p))) / 2) / 4
ab
Substitute old: j - j', new: j_m
Working...
1:  D k k'
*((exp(i x k_m) + exp(-(i x k_m))) / 2
+ (exp(i x k_p) + exp(-(i x k_p))) / 2)
((exp(i y j_m) + exp(-(i y j_m))) / 2
- (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2) / 4
+ D j j' l
*((exp(i y j_m) + exp(-(i y j_m))) / 2
+ (exp(i y*(j + j')) + exp(-(i y*(j + j')))) / 2)
((exp(i x k_m) + exp(-(i x k_m))) / 2
- (exp(i x k_p) + exp(-(i x k_p))) / 2) / 4
ab
Substitute old: j + j', new: j_p
1:  D k k'
*((exp(i x k_m) + exp(-(i x k_m))) / 2
+ (exp(i x k_p) + exp(-(i x k_p))) / 2)
((exp(i y j_m) + exp(-(i y j_m))) / 2
- (exp(i y j_p) + exp(-(i y j_p))) / 2) / 4
+ D j j' l
*((exp(i y j_m) + exp(-(i y j_m))) / 2
+ (exp(i y j_p) + exp(-(i y j_p))) / 2)
((exp(i x k_m) + exp(-(i x k_m))) / 2
- (exp(i x k_p) + exp(-(i x k_p))) / 2) / 4

The next step is to bite the bullet and expand it:
ax
Working...
1:  D k k' exp(i x k_m) exp(i y j_m) / 16
+ D k k' exp(i x k_m) exp(-(i y j_m)) / 16
+ D k k' exp(-(i x k_m)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_m)) exp(-(i y j_m)) / 16
- D k k' exp(i x k_m) exp(i y j_p) / 16
- D k k' exp(i x k_m) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_m)) exp(i y j_p) / 16
- D k k' exp(-(i x k_m)) exp(-(i y j_p)) / 16
+ D k k' exp(i x k_p) exp(i y j_m) / 16
+ D k k' exp(i x k_p) exp(-(i y j_m)) / 16
+ D k k' exp(-(i x k_p)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_p)) exp(-(i y j_m)) / 16
- D k k' exp(i x k_p) exp(i y j_p) / 16
- D k k' exp(i x k_p) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
+ D j j' l exp(i y j_m) exp(i x k_m) / 16
+ D j j' l exp(i y j_m) exp(-(i x k_m)) / 16
+ D j j' l exp(-(i y j_m)) exp(i x k_m) / 16
+ D j j' l exp(-(i y j_m)) exp(-(i x k_m)) / 16
- D j j' l exp(i y j_m) exp(i x k_p) / 16
- D j j' l exp(i y j_m) exp(-(i x k_p)) / 16
- D j j' l exp(-(i y j_m)) exp(i x k_p) / 16
- D j j' l exp(-(i y j_m)) exp(-(i x k_p)) / 16
+ D j j' l exp(i y j_p) exp(i x k_m) / 16
+ D j j' l exp(i y j_p) exp(-(i x k_m)) / 16
+ D j j' l exp(-(i y j_p)) exp(i x k_m) / 16
+ D j j' l exp(-(i y j_p)) exp(-(i x k_m)) / 16
- D j j' l exp(i y j_p) exp(i x k_p) / 16
- D j j' l exp(i y j_p) exp(-(i x k_p)) / 16
- D j j' l exp(-(i y j_p)) exp(i x k_p) / 16
- D j j' l exp(-(i y j_p)) exp(-(i x k_p)) / 16

Now the tricky part is to collect it all into a Fourier transform.

We already have here terms that are proportional to

and these are just fine. But what to do about terms that have a plus in front of one or both is?

• Remember that we're really working with integrals, i.e., we're interested in integrals such as:

We can always change the variables in this integral from x to x' = -x with the following result:

At this stage a lot depends on function D and on its analytical extension into . If D is such that for , then the above integral can be simply translated into without change:

and so we have arrived at

Another way to see the same is to replace and observe that since k+ is an integer .

• So we are justified to carry out the following substitutions:
ab
Substitute old: exp(i x k_p), new: exp(-i x k_p)
Working:
1:  D k k' exp(i x k_m) exp(i y j_m) / 16
+ D k k' exp(i x k_m) exp(-(i y j_m)) / 16
+ D k k' exp(-(i x k_m)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_m)) exp(-(i y j_m)) / 16
- D k k' exp(i x k_m) exp(i y j_p) / 16
- D k k' exp(i x k_m) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_m)) exp(i y j_p) / 16
- D k k' exp(-(i x k_m)) exp(-(i y j_p)) / 16
+ D k k' exp(-(i x k_p)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_p)) exp(-(i y j_m)) / 16
+ D k k' exp(-(i x k_p)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_p)) exp(-(i y j_m)) / 16
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
+ D j j' l exp(i y j_m) exp(i x k_m) / 16
+ D j j' l exp(i y j_m) exp(-(i x k_m)) / 16
+ D j j' l exp(-(i y j_m)) exp(i x k_m) / 16
+ D j j' l exp(-(i y j_m)) exp(-(i x k_m)) / 16
- 1:8 D j j' l exp(i y j_m) exp(-(i x k_p))
- 1:8 D j j' l exp(-(i y j_m)) exp(-(i x k_p))
+ D j j' l exp(i y j_p) exp(i x k_m) / 16
+ D j j' l exp(i y j_p) exp(-(i x k_m)) / 16
+ D j j' l exp(-(i y j_p)) exp(i x k_m) / 16
+ D j j' l exp(-(i y j_p)) exp(-(i x k_m)) / 16
- 1:8 D j j' l exp(i y j_p) exp(-(i x k_p))
- 1:8 D j j' l exp(-(i y j_p)) exp(-(i x k_p))
ab
Substitute old: exp(i x k_m), new: exp(-i x k_m)
Working...
1:  D k k' exp(-(i x k_m)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_m)) exp(-(i y j_m)) / 16
+ D k k' exp(-(i x k_m)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_m)) exp(-(i y j_m)) / 16
- D k k' exp(-(i x k_m)) exp(i y j_p) / 16
- D k k' exp(-(i x k_m)) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_m)) exp(i y j_p) / 16
- D k k' exp(-(i x k_m)) exp(-(i y j_p)) / 16
+ D k k' exp(-(i x k_p)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_p)) exp(-(i y j_m)) / 16
+ D k k' exp(-(i x k_p)) exp(i y j_m) / 16
+ D k k' exp(-(i x k_p)) exp(-(i y j_m)) / 16
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
+ 1:8 D j j' l exp(i y j_m) exp(-(i x k_m))
+ 1:8 D j j' l exp(-(i y j_m)) exp(-(i x k_m))
- 1:8 D j j' l exp(i y j_m) exp(-(i x k_p))
- 1:8 D j j' l exp(-(i y j_m)) exp(-(i x k_p))
+ 1:8 D j j' l exp(i y j_p) exp(-(i x k_m))
+ 1:8 D j j' l exp(-(i y j_p)) exp(-(i x k_m))
- 1:8 D j j' l exp(i y j_p) exp(-(i x k_p))
- 1:8 D j j' l exp(-(i y j_p)) exp(-(i x k_p))
ab
Substitute old: exp(i y j_m), new: exp(-i y j_m)
Working...
1:  1:4 D k k' exp(-(i x k_m)) exp(-(i y j_m))
- D k k' exp(-(i x k_m)) exp(i y j_p) / 16
- D k k' exp(-(i x k_m)) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_m)) exp(i y j_p) / 16
- D k k' exp(-(i x k_m)) exp(-(i y j_p)) / 16
+ 1:4 D k k' exp(-(i x k_p)) exp(-(i y j_m))
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
- D k k' exp(-(i x k_p)) exp(i y j_p) / 16
- D k k' exp(-(i x k_p)) exp(-(i y j_p)) / 16
+ 1:4 D j j' l exp(-(i y j_m)) exp(-(i x k_m))
- 1:4 D j j' l exp(-(i y j_m)) exp(-(i x k_p))
+ 1:8 D j j' l exp(i y j_p) exp(-(i x k_m))
+ 1:8 D j j' l exp(-(i y j_p)) exp(-(i x k_m))
- 1:8 D j j' l exp(i y j_p) exp(-(i x k_p))
- 1:8 D j j' l exp(-(i y j_p)) exp(-(i x k_p))
ab
Substitute old: exp(i y j_p), new: exp(-i y j_p)
Working...
1:  1:4 D k k' exp(-(i x k_m)) exp(-(i y j_m))
- 1:4 D k k' exp(-(i x k_m)) exp(-(i y j_p))
+ 1:4 D k k' exp(-(i x k_p)) exp(-(i y j_m))
- 1:4 D k k' exp(-(i x k_p)) exp(-(i y j_p))
+ 1:4 D j j' l exp(-(i y j_m)) exp(-(i x k_m))
- 1:4 D j j' l exp(-(i y j_m)) exp(-(i x k_p))
+ 1:4 D j j' l exp(-(i y j_p)) exp(-(i x k_m))
- 1:4 D j j' l exp(-(i y j_p)) exp(-(i x k_p))

• This we can already collect with ease. Let us first introduce

 (3.45)

Then the last expression returned by Calc becomes:
 (3.46)

• And so our equation of motion becomes now

• We can sweep the complexity of this equation under the carpet and introduce a multidimensional matrix Fkjk'j' such that:

 (3.47)

What exactly goes into Fkjk'j' can be easily inferred from equation (3.49).
• But this still doesn't look like a simple matrix equation. We can fix that by doing two things:
1.
Truncating the infinite series in equation (3.50) at some finite k'=nx and j'=ny.
2.
Collapsing dual indices k and j into   a single index l:
 l = (j - 1) nx + k (3.48) j = (3.49) k = (3.50)

At the end of this procedure our diffusion equation transformed to the power space looks as follows:

 (3.51)

Observe a very important difference compared to the Jacobi iteration method. When using spectral method the discretisation does not take place in physical space and/or time. The solution whether exact or approximate is perfectly continuous in space and time. On the other hand even an exact solution is always discretised in the power space, i.e., in k and j. It is discretised, but not finite. We have to truncate the power series in order to obtain a numerical solution. This is where the approximation comes in.

Next: The Eigenvalue Problem Up: Spectral Methods Previous: Spectral Methods
Zdzislaw Meglicki
2001-02-26