{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fourier interpolation using FFT"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format='svg'\n",
"\n",
"from numpy import pi,zeros,arange,exp,sin,abs,linspace,dot,outer,real,mod,log\n",
"from numpy import sqrt,sum\n",
"from numpy.fft import fft\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function computes the trigonometric interpolant and plots it."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Wave numbers are arranged as k=[0,1,...,N/2,-N/2+1,-N/2,...,-1]\n",
"def fourier_interp(N,f,ne=500,fig=True):\n",
" if mod(N,2) != 0:\n",
" print(\"N must be even\")\n",
" return\n",
" h = 2*pi/N; x = h*arange(0,N);\n",
" v = f(x);\n",
" v_hat = fft(v)\n",
" k = zeros(N)\n",
" n = int(N/2)\n",
" k[0:n+1] = arange(0,n+1)\n",
" k[n+1:] = arange(-n+1,0,1)\n",
"\n",
" xx = linspace(0.0,2*pi,ne)\n",
" vf = real(dot(exp(1j*outer(xx,k)), v_hat)/N)\n",
" ve = f(xx)\n",
"\n",
" # Plot interpolant and exact function\n",
" if fig:\n",
" plt.plot(x,v,'o',xx,vf,xx,ve)\n",
" plt.legend(('Data','Fourier','Exact'))\n",
" plt.title(\"N = \"+str(N))\n",
"\n",
" errori = abs(vf-ve).max()\n",
" error2 = sqrt(h*sum((vf-ve)**2))\n",
" print(\"Error (max,L2) = \",errori,error2)\n",
" return errori, error2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Infintely smooth, periodic function"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 8.01581023779363e-14 4.559011236168057e-13\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f1 = lambda x: exp(sin(x))\n",
"fourier_interp(24,f1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Infinitely smooth, derivative not periodic"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 0.024794983262922784 0.06907606978750452\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"g = lambda x: sin(x/2)\n",
"e1i,e12 = fourier_interp(24,g)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 0.012409450043334486 0.017665433009611636\n",
"Rate (max,L2) = 0.9986090713584319 1.9672568878077048\n"
]
}
],
"source": [
"e2i,e22 = fourier_interp(48,g,fig=False)\n",
"print(\"Rate (max,L2) = \", log(e1i/e2i)/log(2), log(e12/e22)/log(2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Continuous function"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 0.03195085708656353 0.1114394822310562\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def trihat(x):\n",
" f = 0*x\n",
" for i in range(len(x)):\n",
" if x[i] < 0.5*pi or x[i] > 1.5*pi:\n",
" f[i] = 0.0\n",
" elif x[i] >= 0.5*pi and x[i] <= pi:\n",
" f[i] = 2*x[i]/pi - 1\n",
" else:\n",
" f[i] = 3 - 2*x[i]/pi\n",
" return f\n",
" \n",
"e1i,e12 = fourier_interp(24,trihat)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 0.015806700636915583 0.028104406730080144\n",
"Rate (max,L2) = 1.0153183695992953 1.9873921942542216\n"
]
}
],
"source": [
"e2i,e22 = fourier_interp(48,trihat,fig=False)\n",
"print(\"Rate (max,L2) = \", log(e1i/e2i)/log(2), log(e12/e22)/log(2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Discontinuous function"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 0.9915112319641801 2.0366611449817986\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f2 = lambda x: (abs(x-pi) < 0.5*pi)\n",
"e1i,e12 = fourier_interp(24,f2)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error (max,L2) = 0.9828401566916884 1.041560145620768\n",
"Rate (max,L2) = 0.0126723113206373 0.9674598167513287\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"e2i,e22 = fourier_interp(48,f2)\n",
"print(\"Rate (max,L2) = \", log(e1i/e2i)/log(2), log(e12/e22)/log(2))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}