My assignment is to compare analytical and numerical solutions of a homogenous heat equation $u_t=au_{xx}$ with initial condition
$u(x,0)=\phi(x)=\sin(x)+\frac{1}{10}\sin(10x)$
and boundary conditions
$u(0,t)=0$
$u(l,t)=0$.
The domain to solve in is $x\in<0,2\pi>, t\in<0,T>$ (T wasn't given so I choose $1$). The numerical solution didn't pose a problem, but I'm struggling with the analytical. The expected result is a dissipating sine wave, starting with the initial condition of two superimposed waves, where the higher frequency wave dies out rapidly. Expected solution. At first I tried the Fourier transform method, i.e. $$u(x,t)=\frac{1}{2\sqrt{a\pi t}}\int_{-\infty}^{\infty}\phi(\xi)\exp\left(-\frac{(x-\xi)^2}{4at}\right)d\xi$$
$$u(x,t)=\frac{1}{2\sqrt{a\pi t}}\int_{-\infty}^{\infty}\left(\sin(\xi)+\frac{1}{10}\sin(10\xi)\right)\exp\left(-\frac{(x-\xi)^2}{4at}\right)d\xi$$
To evaluate the integral I used the integral command in Matlab. Aside from the evaluation being very slow (which is to be expected) I got generally correct results, but random spikes appear in the solution. My only guess is that it has to do with the internal workings of the integral command.
Fourier transformation solution
My next try was Fourier series. This resulted in integrals simple enough, that I could solve by hand and avoid the numerical integration.
$$u(x,t)=\sum_{k=1}^{\infty}\frac{2}{l}\int_0^l\phi(\xi)\sin\left(\frac{k\pi\xi}{l}\right)d\xi \exp\left(-\left(\frac{k\pi}{l}\right)^2t\right)\sin\left(\frac{k\pi x}{l}\right)$$
$$u(x,t)=\sum_{k=1}^{\infty}\frac{2}{l}\int_0^l\left(\sin(\xi)+\frac{1}{10}\sin(10\xi)\right)\sin\left(\frac{k\pi\xi}{l}\right)d\xi \exp\left(-\left(\frac{k\pi}{l}\right)^2t\right)\sin\left(\frac{k\pi x}{l}\right)$$
$$u(x,t)=\sum_{k=1}^{\infty}\frac{2}{l}\left(I+II\right)\exp\left(-\left(\frac{k\pi}{l}\right)^2t\right)\sin\left(\frac{k\pi x}{l}\right)$$
where
$$I=\frac{\sin\left(1-\frac{k\pi}{l}\right)}{2\left(1-\frac{k\pi}{l}\right)}-\frac{\sin\left(1+\frac{k\pi}{l}\right)}{2\left(1+\frac{k\pi}{l}\right)}\hspace{1cm}if:\:k\neq\frac l \pi$$
$$II=\frac {1}{10} \left(\frac{\sin\left(10-\frac{k\pi}{l}\right)}{2\left(10-\frac{k\pi}{l}\right)}-\frac{\sin\left(10+\frac{k\pi}{l}\right)}{2\left(10+\frac{k\pi}{l}\right)}\right)\hspace{1cm}if:\:k\neq\frac{10l}{\pi}$$
$$I=\frac 1 2 l-\frac 1 4 \sin{2l}\hspace{1cm}if:\:k=\frac l \pi$$
$$II=\frac {1}{10}\left(\frac 1 2 l-\frac{1}{40}\sin{20l}\right)\hspace{1cm}if:\:k=\frac{10l}{\pi}$$
Unfortunately, the results are way off. Fourier series solution Not really sure what's the problem here, the coefficents seem correct to me.
At last, I tried to use FFT. We didn't learn this method, so I mainly followed a YT video. It looked rather simple, I know the solution of the transformed equation ($y$ is the angular frequency):
$$\hat u(y,t)=\hat \phi(y)\exp(-ay^2t)$$
I know I need to transform the initial condition, what I don't understand is what should the angular frequency values be. I understand the reason behind $2\pi$ as $\omega=2\pi\xi$ ($\xi$ is ordinary frequency) but not much else.
No picture to show here as the inverse transformation always returns complex numbers.
Any advice here would be greatly appreciated. Not sure I can do much about the Fourier transform solution, but I had high hopes for the series expansion. Also the FFT seemed simple enough. But I must be missing a key idea somewhere.
Here's the Matlab code for all attempts at solution.
Notes: I showed only solutions for $a=1$ but the same issues arise for $a=10, 100$
The reason for such tight time discretization is the stability of numerical solutions. I know I could increase the time step, but with the same mesh I can compare the solutions easily.
I also tried Matlab's pdepe command, which gave correct solution, but as I haven't provided the solution, I can't present it to my teacher.
Best Answer
Using Fourier series
Since we have $u_{xx}$ and boundary conditions $u(0,t)=0,\ u(2\pi,t)=0$ we expand $u(x,t)$ in the basis $\{ \sin\frac{kx}{2} : k=1,2,3,\ldots\}$: $$u(x,t) = \sum_{k=1}^{\infty} \hat{u}_k(t) \sin\frac{kx}{2}.$$
Inserting this into the PDE gives $$ \sum_{k=1}^{\infty} a\hat{u}_k'(t) \sin\frac{kx}{2} = -\sum_{k=1}^{\infty} \hat{u}_k(t) \left(\frac{k}{2}\right)^2 \sin\frac{kx}{2} $$ from which we conclude that for each $k$ and all $t>0,$ $$ a\hat{u}_k'(t) = -\left(\frac{k}{2}\right)^2 \hat{u}_k(t) . $$ Thus, $$ \hat{u}_k(t) = \hat{u}_k(0) e^{-k^2t/(4a)} $$ and $$ u(x,t) = \sum_{k=1}^{\infty} \hat{u}_k(0) e^{-k^2t/(4a)} \sin\frac{kx}{2} . $$
We also have an initial condition $u(x,0)=\sin(x)+\frac{1}{10}\sin(10x),$ which implies that $\hat{u}_2(0)=1$ and $\hat{u}_{20}(0)=\frac{1}{10}$ while $\hat{u}_k(0)=0$ for all other values of $k$. Thus we end up with $$ u(x,t) %= \sum_{k\in\{2,20\}} \hat{u}_k(0) e^{-k^2t/(4a)} \sin\frac{kx}{2} = e^{-t/a} \sin x + \frac{1}{10} e^{-t/(100a)} \sin 100x . $$