I'm trying to solve a one-dimensional heat equation (PDE) with the Fourier transform numerically, in the way it was done here. The equation:
is subject to the initial condition:
where U(x,t) is temperature, x is space, a is heat conductivity, and t is time.
I want to solve this equation using fast Fourier transform (FFT).
The key observation here is concerning the derivatives:
where k=2 pi/L[-N/2,N/2] is a spatial frequency or wave number. So, u(k,t) is a vector of Fourier coefficients and k square is a vector of frequency, so that
gives n decoupled ODEs, one for each these kj.
The respective code is:
n = 1000;(*Number of discretization points*)
L = 100; (*Length of domain**)
T = 20; (*Time Integration*)
a = 1; (*thermal diffusivity constant*)
k = (2 Pi)/L Table[i, {i, -n/2, n/2 - 1}];(*define discrete wave number*)
kt = Fourier[k];(*Re-order FFT wave number*)
ic1 = Table[If[400 < i < 600, 1, 0], {i, Length[k]}]; (*initial condition*)
ict = Fourier[ic1];(*Fourier Transform of initial condition*)
ic = Table[Subscript[u, i][0] == ict[[i]], {i,Length[k]}];(*vector of initial condition in Fourier transform domain*)
vars = Table[Subscript[u, i][t], {i, Length[k]}]; (*vector of variables*)
eqns = Table[{Subscript[u, i]'[t] == -a (kt[[i]])^2 Subscript[u, i][t]}, {i,Length[k]}];(*model of ODEs system*)
eqn = Join[eqns, ic];
sol = NDSolve[eqn,vars, {t, 0, T}];(*Simulate in Fourier frequency domain*)
This error message appears:
General::munfl: (-1.4792510^-23+4.7229210^-303> I)+(1.5440710^-23-4.7229110^-303 I) is too small to represent as a> normalized machine number; precision may be lost. General::stop:
Further output of General::munfl will be suppressed during this > calculation.
How can I fix this?
Next, I intend to plot the solution by applying the inverse Fourier transform (IFFT).
s = Table[ Table[Evaluate[Subscript[u, j][t] /. sol], {t, 0, T}], {j,Length[k]}];(*select the solution of temperature in frequency domain*)
sin = InverseFourier[s];(*IFFT to return to spatial domain*)
fin = Table[{i, j, sin[[j, i, 1]]}, {i, Length[sin[[1]]]}, {j,Length[sin]}];
Table[ListPlot3D[fin[[i]]], {i, Length[fin]}];
The plot of the solution I'm looking for should be something like: