A Diffusion Problem

• Consider the following slightly more complex diffusion equation:
  

  $$\frac{\partial}{\partial t} T(x, y, t) = \boldsymbol{\nabla}\cdot D(x, y) \boldsymbol{\nabla} T(x, y, t)$$ (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 l_x and l_y.
• We can always   change our units of length and time so that the beam cross section becomes a square with a side length of $\pi$ - this will lead to certain FFT conveniences. The transformation of coordinates is:
  

  $$\pi x / l_x$$ x' = (3.25)

  $$\pi y / l_y$$ y' = (3.26)

  
  

• This implies
  

  $$\frac{\partial}{\partial x} \frac{\pi}{l_x}\frac{\partial}{\partial x'}$$ = (3.27)

  $$\frac{\partial}{\partial y} \frac{\pi}{l_y}\frac{\partial}{\partial y'}$$ = (3.28)

  
  

• Substitute into equation (3.24) to get:
  

  $$\frac{\partial}{\partial t} T \left(\frac{\pi}{l_x}\right)^2 \left( \frac{\partial}{\partial x'... ...ial y'} \left(\frac{l_x}{l_y}\right)^2 D \frac{\partial}{\partial y'} \right) T$$ =

  $$\left(\frac{\pi}{l_x}\right)^2 \left( \frac{\partial}{\partial x'... ...} + \frac{\partial}{\partial y'} l_r^2 D \frac{\partial}{\partial y'} \right) T$$ = (3.29)

  
  
  where l_r = l_x/l_y.
• The Diffusion coefficient   has units of $\hbox{m}^2\,\hbox{s}^{-1}$. 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, D_r and D itself will become D'. In summary D = D_r D', where D_r is measured in $\hbox{m}^2$ and:
  

  $$\frac{\partial}{\partial t} T = \frac{\pi^2 D_r}{l_x^2} \lef... ...rtial y'} l_r^2 D' \frac{\partial}{\partial y'} \right) T$$ (3.30)

  
  

• The factor $\pi^2 D_r / l_x^2$ can be now absorbed into a new unit of time (observe that $\hbox{m}^2$ of D_r cancels with l_x² so that the new unit of time still has the corect dimension of some fraction of seconds):
  

  $$\frac{l_x^2}{D_r\pi^2} dt = dt'$$ (3.31)

  
  
  so that
  

  $$\frac{\partial}{\partial t'} T = \left( \frac{\partial}{... ...rtial y'} l_r^2 D' \frac{\partial}{\partial y'} \right) T$$ (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 $\pi$ 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) = T₀(x, y) (3.33)

  
  

• Assume the following fixed boundary conditions:  
  

  $$T(0, y, t) = T(\pi, y, t) = T(x, 0, t) = T(x, \pi, t) = 0$$ (3.34)

  
  

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

  $$T(x, y, t) = \frac{2}{\pi} \sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j'}(t)\sin(k'x)\sin(j'y)$$ (3.35)

  
  
  because for x=0 or $x=\pi$ $\sin(k'x) = 0$ and for y=0 or $y=\pi$ $\sin(j'y) = 0$ too. The factor $\frac{2}{\pi}$ is there for future convenience. You can always shove it into a_k'j' if you don't like it.
• Given T(x, y, 0) = T₀(x, y) it is easy to evaluate a_k'j'(0). First observe that
  

  $$\int_0^\pi \sin(n x) \sin(m x) \,\textrm{d} x = \frac{\pi}{2}\delta_{nm}$$ (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 $\pi/2$.

  We make use of the above as follows:
  

  $$\frac{2}{\pi}\int_0^\pi\int_0^\pi \sin(kx)\sin(jy)T_0(x,y)\,\textrm{d}x\,\textrm{d}y$$

  $$\quad = \frac{2}{\pi}\int_0^\pi\int_0^\pi \sin(kx)\sin(jy)\frac{2... ...infty\sum_{j'=1}^\infty a_{k'j'}(0)\sin(k'x)\sin(j'y)\,\textrm{d}x\,\textrm{d}y$$

  $$\quad = \frac{4}{\pi^2}\sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'... ...nt_0^\pi\int_0^\pi\sin(kx)\sin(k'x)\sin(jy)\sin(j'y) \,\textrm{d}x\,\textrm{d}y$$

  $$\quad = \frac{4}{\pi^2}\sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j'}(0) \frac{\pi}{2}\delta_{kk'}\frac{\pi}{2}\delta_{jj'}$$

  $$\quad = a_{kj}(0)$$ (3.37)

  
  
  Now you can also see why we have defined our solution with that $2/\pi$ 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
  

  $$\int_0^{2\pi} \cos(m x) \sin(n x)\,\textrm{d} x = \frac{1}{n... ...1\right)} + \frac{1}{n\left(1 - \frac{m^2}{n^2}\right)} = 0$$ (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:
  

  $$\int_0^{2\pi} \cos(n x) \sin(n x)\,\textrm{d} x = 0$$ (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:
  

  $$\int_0^{2\pi} \sin(m x) \sin(n x)\,\textrm{d} x = 0$$ (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 $\pi$:
  

  $$\int_0^{2\pi} \sin(n x) \sin(n x)\,\textrm{d} x = \pi$$ (3.41)

  
  

• We can now extend T₀(x, y) analytically to a $[0, 2\pi]\times[0, 2\pi]$ region by keeping the same coefficients, but letting x and y run between 0 and $2\pi$ each.  
• Now consider the following integral:
  
  $$\begin{eqnarray*}&&\frac{-1}{2\pi} \int_0^{2\pi}\int_0^{2\pi}e^{-ikx}e^{-ijy}T_... ...y)\right) \sin(k'x)\sin(j'y) \,\textrm{d} x\,\textrm{d} y \\ \end{eqnarray*}$$
  
  This will produce a variety of terms. First there will be terms containing $\cos(kx)\sin(k'x)$ or $\cos(jy)\sin(j'y)$ and they will all vanish on integration over $[0, 2\pi]$, 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 $i^2\sin(kx)\sin(jy)$:
  

  $$\frac{-1}{2\pi} \int_0^{2\pi}\int_0^{2\pi}e^{-ikx}e^{-ijy}T_0(x,y) \,\textrm{d} x\,\textrm{d} y$$

  $$\quad = \frac{-1}{\pi^2} \sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{... ...nt_0^{2\pi} - \sin(kx) \sin(jy) \sin(k'x)\sin(j'y) \,\textrm{d} x\,\textrm{d} y$$

  $$\quad = \frac{1}{\pi^2} \sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j'}(0) \int_0^{2\pi}\int_0^{2\pi} \pi\delta_{kk'}\pi\delta_{jj'}$$

  $$\quad = a_{kj}(0)$$ (3.42)

  
  
   
• Let us now plug T(x, y, t) as given by quation (3.37) into our diffusion equation (3.34):
  
  $$\begin{eqnarray*}&&\frac{2}{\pi}\sum_{k'=1}^\infty\sum_{j'=1}^\infty \sin(k' x)... ...}^\infty a_{k'j'}(t) \hat{\boldsymbol{D}} \sin(k' x)\sin(j' y) \end{eqnarray*}$$
  
  where $\hat{\boldsymbol{D}} = \partial_x D \partial_x + \partial_y l_r^2 D \partial_y$.
• The next step is to multiply both sides by $\frac{2}{\pi}\sin(kx)\sin(jy)$ and integrate over $[0,\pi]$. 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:
  
  $$\begin{eqnarray*}&& \int_0^\pi\int_0^\pi \frac{2}{\pi}\sin(kx)\sin(jy) \frac{2... ...{jj'} \\ &&\quad = \frac{\textrm{d} a_{kj}(t)}{\textrm{d} t} \end{eqnarray*}$$
  
  The right hand side is a little more complicated because of the presence of differential operators and an unknown function D(x,y):
  
  $$\begin{eqnarray*}&& \int_0^\pi\int_0^\pi \frac{2}{\pi}\sin(kx)\sin(jy) \frac{2... ...boldsymbol{D}} \sin(k' x)\sin(j' y) \,\textrm{d}x\,\textrm{d}y \end{eqnarray*}$$
  
  

• Before going any further let us evaluate
  
  $$\hat{\boldsymbol{D}}\sin(k'x)\sin(j'y)$$
  
  
  first. It is easy to see that this is:
  
  $$\frac{\partial}{\partial x} \left( D k'\cos(k'x)\sin(j'y) ... ...l}{\partial y} \left( l_r^2 D j'\sin(k'x)\cos(j'y) \right)$$
  
  

• Now plug this into the integral to get this very lengthy expression:
  
  $$\begin{eqnarray*}&& \frac{4}{\pi^2}\sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j... ...sin(k'x)\cos(j'y) \right) \right) \,\textrm{d}x\,\textrm{d}y \end{eqnarray*}$$
  
  
• We can't differentiate D because we don't know what it is, but we can integrate by parts and transfer $\partial/\partial x$ and $\partial/\partial y$ onto $\sin(kx)\sin(jy)$. The result is:
  

  $$\frac{\textrm{d} a_{kj}(t)}{\textrm{d} t} =$$

  $$\quad - \frac{4}{\pi^2}\sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j'}(t) \int_0^\pi\int_0^\pi \big( kk'D \cos(kx)\cos(k'x)\sin(jy)\sin(j'y)$$

  $$\qquad\qquad + \> l_r^2 jj'D \cos(jy)\cos(j'y)\sin(kx)\sin(k'x) \big) \,\textrm{d}x\,\textrm{d}y$$ (3.43)

  
  

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

  $$\frac{\textrm{d} a_{kj}(t)}{\textrm{d} t} =$$

  $$\quad - \frac{1}{\pi^2}\sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j'}(t) \int_0^{2\pi}\int_0^{2\pi} \big( kk'D \cos(kx)\cos(k'x)\sin(jy)\sin(j'y)$$

  $$\qquad\qquad + \> l_r^2 jj'D \cos(jy)\cos(j'y)\sin(kx)\sin(k'x) \big) \,\textrm{d}x\,\textrm{d}y$$ (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:
    
    $$\begin{eqnarray*}\cos(k x)\cos(k' x) &=& \frac{1}{2}\left( \cos\left((k - k') x... ...left((k - k') x\right) - \cos\left((k + k'\right) x) \right) \end{eqnarray*}$$
    
     

    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:
    
    $$\cos(kx) = \frac{1}{2}\left(e^{i k x} + e^{-i k x}\right)$$
    
    

    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
    
    $$\begin{eqnarray*}k + k' &=& k^+ \\ k - k' &=& k^- \\ j + j' &=& j^+ \\ j - j' &=& j^- \\ \end{eqnarray*}$$
    
      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
    

    $$e^{-i k^\pm x}e^{-i j^\pm y}$$
    
    
    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:
    
    $$\int_0^{2\pi} D(x, y) e^{i k^+ x} \,\textrm{d} x$$
    
    
    We can always change the variables in this integral from x to x' = -x with the following result:
    
    $$- \int_0^{-2\pi} D(x', y) e^{-i k^+ x'}\,\textrm{d} x' = \int_{-2\pi}^0 D(x', y) e^{-i k^+ x'}\,\textrm{d} x'$$
    
    
    At this stage a lot depends on function D and on its analytical extension into $[-2\pi, 0]$. If D is such that $D(-2\pi + x) = D(x)$ for $x\in[0, 2\pi]$, then the above integral can be simply translated into $[0, 2\pi]$ without change:
    
    $$\int_{-2\pi}^0 D(x', y) e^{-i k^+ x'}\,\textrm{d} x' = \int_0^{2\pi} D(x, y) e^{-i k^+ x}\,\textrm{d} x$$
    
    
    and so we have arrived at
    
    $$\int_0^{2\pi} D(x, y) e^{i k^+ x}\,\textrm{d} x = \int_0^{2\pi} D(x, y) e^{-i k^+ x}\,\textrm{d} x$$
    
    
    Another way to see the same is to replace $x' = 2\pi - x$ and observe that since k⁺ is an integer $e^{-i k 2\pi} = 1$.

  • 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
    

    $$\tilde{D}_{kj} = D e^{-i x k} e^{-i y j}$$ (3.45)

    
    
    Then the last expression returned by Calc becomes:
    

    $$\frac{1}{4} k k' \left( \tilde{D}_{k^-j^-} - \tilde{D}_{k^-j^+} + \tilde{D}_{k^+j^-} - \tilde{D}_{k^+j^+} \right)$$

    $$\quad + \frac{1}{4}l^2_r j j' \left( \tilde{D}_{k^-j^-} - \tilde{D}_{k^+j^-} + \tilde{D}_{k^-j^+} - \tilde{D}_{k^+j^+} \right)$$

    $$= \frac{1}{4} \bigg( \left( k k' + l^2_r j j' \right) \left( \tilde{D}_{k^-j^-} - \tilde{D}_{k^+j^+} \right)$$

    $$\quad + \left( k k' - l^2_r j j' \right) \left( \tilde{D}_{k^+j^-} - \tilde{D}_{k^-j^+} \right) \bigg)$$ (3.46)

    
    

• And so our equation of motion becomes now  
  
   

  $$\frac{\textrm{d} a_{kj}(t)}{\textrm{d} t} = - \frac{1}{4\pi^2}\sum_{k'=1}^\infty\sum_{j'=1}^\infty a_{k'j'}(t)$$

  $$\times \bigg( \left(k k' + l_r^2 j j'\right) \left( \int_0^{2\pi}... ... \int_0^{2\pi}\int_0^{2\pi}\tilde{D}_{k^+j^+}\,\textrm{d}x\,\textrm{d}y \right)$$

  $$+ \left(k k' - l_r^2 j j'\right) \left( \int_0^{2\pi}\int_0^{2\pi... ...^{2\pi}\int_0^{2\pi}\tilde{D}_{k^-j^+}\,\textrm{d}x\,\textrm{d}y \right) \bigg)$$

  
  

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

  $$\frac{\textrm{d} a_{kj}(t)}{\textrm{d} t} = -\sum_{k'j'}F_{kjk'j'}a_{k'j'}(t)$$ (3.47)

  
  
  What exactly goes into F_kjk'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'=n_x and j'=n_y.

  2.

  Collapsing dual indices k and j into   a single index l:
  

   

  l = (j - 1) n_x + k (3.48)

  $$\frac{l - 1}{n_x} + 1$$ j = (3.49)

  $$\hbox{mod}(l - 1, n_x) + 1$$ k = (3.50)

  
  

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

  $$\frac{\textrm{d} a_l(t)}{\textrm{d} t} = -\sum_{l'=1}^{n_x n_y}F_{ll'}a_{l'}(t)$$ (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.