Diffusion through Fourier modes.
Watch short wavelengths disappear first.
The physical idea.
For ∂c/∂t = D∇²c, each Fourier amplitude decays as ĉ(k,t) = ĉ(k,0) exp(−Dk²t). Here c(x,0) = 1 + 0.4 cos(x) + 0.2 cos(4x), D = 1 and the periodic domain is [0,2π).
This example is based on my teaching notebook, Fourier_Julia.ipynb. I wrote the notebook to explain how Fourier methods solve the diffusion equation. Here, two Fourier modes show why short-wavelength variations disappear faster.
Precomputed NumPy results from the downloadable code; the browser scrubs saved states. Line: concentration versus x.
Download measured diagnostics (CSV) ↗Model and numerical choices.
This is exact propagation of the represented Fourier modes. There is no explicit time-step stability limit. The fourth harmonic decays sixteen times faster than the first. Concentration and coordinates are dimensionless.
Questions explored in this example.
Change the fourth harmonic to an eighth harmonic. Predict its decay rate before running the code. Verify that mean concentration stays at 1.
My Python calculation, explained in steps.
Both versions implement the same model. Python uses NumPy; Julia uses FFTW for spectral examples. Scripts write CSV snapshots and diagnostics into the output folder.
Download the Python notebook ↓ · Read the code in steps → Download Julia version ↓
Run the notebook cells in order for the explained calculation. The commands and complete script below are an optional command-line route; the script is linked at the end of the notebook lesson.
python -m pip install numpy
python diffusion.py resultsRead the complete NumPy code
"""Diffusion example based on my teaching notebook Fourier_Julia.ipynb.
Purpose: show how Fourier modes decay during diffusion.
Run: python diffusion.py [output_directory]. Requires NumPy.
Periodic domain L=2*pi; D=1. Exact Fourier propagation, not time stepping.
"""
import sys
from pathlib import Path
import numpy as np
def run(output='diffusion_output'):
out=Path(output);out.mkdir(parents=True,exist_ok=True)
n=128;x=2*np.pi*np.arange(n)/n;k=2*np.pi*np.fft.fftfreq(n,d=2*np.pi/n)
c0=1+.4*np.cos(x)+.2*np.cos(4*x);rows=[]
for t in [0,.02,.1,.3,1.]:
c=np.fft.ifft(np.fft.fft(c0)*np.exp(-k*k*t)).real
exact=1+.4*np.exp(-t)*np.cos(x)+.2*np.exp(-16*t)*np.cos(4*x)
rows.append([t,c.mean(),np.max(np.abs(c-exact))])
np.savetxt(out/f'profile_{t:.2f}.csv',np.column_stack([x,c]),delimiter=',')
np.savetxt(out/'diagnostics.csv',rows,delimiter=',',header='time,mean,max_exact_error',comments='')
assert max(r[2] for r in rows)<1e-12
print('Fourier diffusion: max exact error',max(r[2] for r in rows))
if __name__=='__main__':run(*sys.argv[1:2])
julia -e 'using Pkg; Pkg.add("FFTW")'
julia diffusion.jl resultsRead the complete Julia code
# Diffusion example based on my teaching notebook Fourier_Julia.ipynb.
# Purpose: show how Fourier modes decay during diffusion.
# julia diffusion.jl [output_directory]; requires FFTW.
using FFTW, Statistics, DelimitedFiles, Printf
function run(output="diffusion_output")
mkpath(output);n=128;x=2π .* (0:n-1) ./ n;k=2π .* fftfreq(n,n/(2π))
c0=1 .+ .4 .* cos.(x) .+ .2 .* cos.(4 .* x);rows=Vector{Float64}[]
for t in [0.,.02,.1,.3,1.]
c=real.(ifft(fft(c0) .* exp.(-k .^ 2 .* t)))
exact=1 .+ .4*exp(-t) .* cos.(x) .+ .2*exp(-16*t) .* cos.(4 .* x)
push!(rows,[t,mean(c),maximum(abs.(c .- exact))])
writedlm(joinpath(output,@sprintf("profile_%.2f.csv",t)),hcat(x,c),',')
end
open(joinpath(output,"diagnostics.csv"),"w") do io
println(io,"time,mean,max_exact_error");writedlm(io,permutedims(hcat(rows...)),',')
end
@assert maximum(r[3] for r in rows)<1e-12
println("Fourier diffusion: max exact error ",maximum(r[3] for r in rows))
end
if abspath(PROGRAM_FILE) == @__FILE__
run(ARGS...)
end
About these calculations.
I use these examples to explain numerical methods. I provide the parameters, boundary conditions and numerical checks with the results. These are teaching calculations, not predictions calibrated for a particular material.
← All worked examplesHow I develop my codes.
I start with a small calculation, check its result and then add the next part. These lessons follow the same sequence, with explanations and checks at each stage. Open the worked Python, Julia and C lessons →