#!/usr/bin/env python # coding: utf-8 # # Fourier spectral methods in Matlab (and Python) # Developed by Randy LeVeque for a course on Approximation Theory and Spectral Methods at the University of Washington. # # See for more IPython Notebook examples. # These examples are based on material in Nick Trefethen's book Spectral Methods in Matlab. The m-files for this book are available at # In[1]: get_ipython().run_line_magic('load_ext', 'pymatbridge') # ## Program 5 # This example is directly from p5.m found at # In[2]: get_ipython().run_cell_magic('matlab', '', '\n% p5.m - repetition of p4.m via FFT\n% For complex v, delete "real" commands.\n\n% Differentiation of a hat function:\n N = 24; h = 2*pi/N; x = h*(1:N)\';\n v = max(0,1-abs(x-pi)/2); v_hat = fft(v);\n w_hat = 1i*[0:N/2-1 0 -N/2+1:-1]\' .* v_hat;\n w = real(ifft(w_hat)); clf\n subplot(2,2,1), plot(x,v,\'.-\',\'markersize\',13)\n axis([0 2*pi -.5 1.5]), grid on, title(\'function\')\n subplot(2,2,2), plot(x,w,\'.-\',\'markersize\',13)\n axis([0 2*pi -1 1]), grid on, title(\'spectral derivative\')\n\n% Differentiation of exp(sin(x)):\n v = exp(sin(x)); vprime = cos(x).*v;\n v_hat = fft(v);\n w_hat = 1i*[0:N/2-1 0 -N/2+1:-1]\' .* v_hat;\n w = real(ifft(w_hat));\n subplot(2,2,3), plot(x,v,\'.-\',\'markersize\',13)\n axis([0 2*pi 0 3]), grid on\n subplot(2,2,4), plot(x,w,\'.-\',\'markersize\',13)\n axis([0 2*pi -2 2]), grid on\n error = norm(w-vprime,inf);\n text(2.2,1.4,[\'max error = \' num2str(error)])\n') # ## Illustration of spectral differentiation # To make this a bit clearer, first illustrate how to compute the second derivative of periodic function. # Start with $$u = \exp(\cos(x)),$$ and check that the numerical approximation agrees well with $$u''(x) = (\sin^2(x) - \cos(x)) \exp(\cos(x)).$$ # # The only tricky thing here is the order of the indices in the wave number vector. # In[3]: get_ipython().run_cell_magic('matlab', '', 'N = 16;\nx = linspace(2*pi/N,2*pi,N);\nik = 1i*[0:N/2 -N/2+1:-1]; % i * wave number vector (matlab ordering)\nik2 = ik.*ik; % multiplication factor for second derivative\n\nu = exp(cos(x));\nu_hat = fft(u);\nv_hat = ik2 .* u_hat;\nv = real(ifft(v_hat)); % imaginary part should be at machine precision level\n\nerror = v - (sin(x).^2 - cos(x)) .* exp(cos(x));\nnorm(error,inf)\n') # ## Illustration of solving a periodic boundary value problem # Now let's solve the boundary value problem # $$u''(x) = f(x)$$ # on $0 \leq x \leq 2\pi$ with periodic boundary conditions and the constraint $\int_0^{2\pi} u(x) dx = 0$. # # Use $f(x) = (\sin^2(x) - \cos(x)) \exp(\cos(x))$ so the solution should be $u(x) = \exp(\cos(x)) + C$, where the constant is chosen so the integral constraint is satisfied. # # We now have to divide by `ik2`, with the complication that 1/0 should be replaced by 0. This results in the $\hat u_0 = 0$, which gives the integral constraint. # In[4]: get_ipython().run_cell_magic('matlab', '', "\nN = 16;\nx = linspace(2*pi/N,2*pi,N);\nf = (sin(x).^2 - cos(x)) .* exp(cos(x));\nf_hat = fft(f);\n\nik = 1i*[0:N/2 -N/2+1:-1]; % i * wave number vector (matlab ordering)\nik2 = ik.*ik; % multiplication factor for second derivative\nii = find(ik ~= 0); % indices where ik is nonzero\nik2inverse = ik2; % initialize zeros in same locations as in ik2\nik2inverse(ii) = 1./ik2(ii); % multiplier factor to solve u'' = f\n\nu_hat = ik2inverse .* f_hat;\nu = real(ifft(u_hat)); % imaginary parts should be roundoff level\n") # Plotting the solution shows that it is a shifted version of $\exp(\cos(x))$: # In[5]: get_ipython().run_cell_magic('matlab', '', "plot(x,u,'b-o')\nhold on\nv = exp(cos(x));\nplot(x,v,'r-o')\n") # If we shift so that one value of $u$ agrees with $v$, then we hope everything will line up: # In[6]: get_ipython().run_cell_magic('matlab', '', 'u2 = u + v(1)-u(1);\nnorm(u2 - v, inf)\n') # ## Python versions: # In[7]: get_ipython().run_line_magic('pylab', 'inline') # We repeat these examples in Python. The codes are essentially identical, with some changes from Matlab to Python notation. # # First illustrate how to compute the second derivative of periodic function. # Start with $$u = \exp(\cos(x)),$$ and check that the numerical approximation agrees well with $$u''(x) = (\sin^2(x) - \cos(x)) \exp(\cos(x))$$ # In[8]: from scipy import fft,ifft N = 16; x = linspace(2*pi/N,2*pi,N) ik = 1j*hstack((range(0,N/2+1), range(-N/2+1,0))); # i * wave number vector (matlab ordering) ik2 = ik*ik; # multiplication factor for second derivative u = exp(cos(x)) u_hat = fft(u) v_hat = ik2 * u_hat v = real(ifft(v_hat)) # imaginary part should be at machine precision level error = v - (sin(x)**2 - cos(x)) * exp(cos(x)) norm(error,inf) # Now let's solve the boundary value problem # $$u''(x) = f(x)$$ # on $0 \leq x \leq 2\pi$ with periodic boundary conditions and the constraint $\int_0^{2\pi} u(x) dx = 0$. # # Use $f(x) = (\sin^2(x) - \cos(x)) \exp(\cos(x))$ so the solution should be $u(x) = \exp(\cos(x)) + C$, where the constant is chosen so the integral constraint is satisfied. # In[9]: N = 16; x = linspace(2*pi/N,2*pi,N) f = (sin(x)**2 - cos(x)) * exp(cos(x)) f_hat = fft(f) ik = 1j*hstack((range(0,N/2+1), range(-N/2+1,0))); # i * wave number vector (matlab ordering) ik2 = ik*ik; # multiplication factor for second derivative ik2inverse = where(ik2 != 0, 1./ik2, 0.) u_hat = ik2inverse * f_hat; u = real(ifft(u_hat)) plot(x,u,'b-o') v = exp(cos(x)); plot(x,v,'r-o') # Again we get good agreement if we shift by the difference at the left-most point: # In[10]: u2 = u + v[0]-u[0] norm(u2 - v, inf) # In[ ]: