gibson:teaching:spring-2014:iam950:hw1
====== Differences ====== This shows you the differences between two versions of the page.
| Both sides previous revision Previous revision Next revision | Previous revision | ||
|
gibson:teaching:spring-2014:iam950:hw1 [2014/02/27 07:32] gibson |
gibson:teaching:spring-2014:iam950:hw1 [2014/02/27 09:39] (current) gibson |
||
|---|---|---|---|
| Line 57: | Line 57: | ||
| with dt = 1/8, 1/4, .... What happens? Why? | with dt = 1/8, 1/4, .... What happens? Why? | ||
| - | **Some help on Matlab's FFT:** Matlab documentation is not very helpful about the ordering or normalization of Fourier coefficients in the FFT. However, a small numerical experiment clarifies: | + | **Some help on Matlab's FFT:** Matlab documentation is not very helpful about the ordering or normalization of Fourier coefficients in the FFT. However, a small numerical experiment clarifies. Let's construct a function with energy in known wavenumbers and observe where nonzero coefficients appear in the numerical FFT. I'll do 8 gridpoints for the periodic domain [0, 2pi]. |
| <code> | <code> | ||
| Line 63: | Line 63: | ||
| x = | x = | ||
| 0 0.7854 1.5708 2.3562 3.1416 3.9270 4.7124 5.4978 | 0 0.7854 1.5708 2.3562 3.1416 3.9270 4.7124 5.4978 | ||
| - | >> u = cos(0*x); | + | |
| - | >> [real(fft(u)); imag(fft(u))] | + | |
| - | ans = | + | |
| - | 8 0 0 0 0 0 0 0 | + | |
| - | 0 0 0 0 0 0 0 0 | + | |
| - | >> u = cos(1*x); | + | |
| - | >> [real(fft(u)); imag(fft(u))] | + | |
| - | ans = | + | |
| - | -0.0000 4.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 4.0000 | + | |
| - | 0 -0.0000 0 -0.0000 0 0.0000 0 0.0000 | + | |
| - | >> u = cos(1*x) + 0.1*sin(1*x); | + | |
| - | >> [real(fft(u)); imag(fft(u))] | + | |
| - | ans = | + | |
| - | -0.0000 4.0000 0.0000 0 0.0000 0 0.0000 4.0000 | + | |
| - | 0 -0.4000 0 -0.0000 0 0.0000 0 0.4000 | + | |
| - | >> | + | |
| - | >> | + | |
| - | >> x = 1/8*(0:7)*2*pi | + | |
| - | x = | + | |
| - | 0 0.7854 1.5708 2.3562 3.1416 3.9270 4.7124 5.4978 | + | |
| >> k = 0; | >> k = 0; | ||
| >> u = cos(k*x) + 0.1*sin(k*x); | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| Line 89: | Line 70: | ||
| 8 0 0 0 0 0 0 0 | 8 0 0 0 0 0 0 0 | ||
| 0 0 0 0 0 0 0 0 | 0 0 0 0 0 0 0 0 | ||
| + | |||
| + | |||
| >> k = 1; | >> k = 1; | ||
| >> u = cos(k*x) + 0.1*sin(k*x); | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| Line 95: | Line 78: | ||
| -0.0000 4.0000 0.0000 0 0.0000 0 0.0000 4.0000 | -0.0000 4.0000 0.0000 0 0.0000 0 0.0000 4.0000 | ||
| 0 -0.4000 0 -0.0000 0 0.0000 0 0.4000 | 0 -0.4000 0 -0.0000 0 0.0000 0 0.4000 | ||
| + | |||
| + | |||
| >> k = 2; | >> k = 2; | ||
| >> u = cos(k*x) + 0.1*sin(k*x); | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| Line 101: | Line 86: | ||
| -0.0000 -0.0000 4.0000 0.0000 0.0000 0.0000 4.0000 -0.0000 | -0.0000 -0.0000 4.0000 0.0000 0.0000 0.0000 4.0000 -0.0000 | ||
| 0 0 -0.4000 0 0 0 0.4000 0 | 0 0 -0.4000 0 0 0 0.4000 0 | ||
| + | |||
| + | |||
| >> k = 3; | >> k = 3; | ||
| >> u = cos(k*x) + 0.1*sin(k*x); | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| Line 107: | Line 94: | ||
| -0.0000 0.0000 -0.0000 4.0000 0.0000 4.0000 -0.0000 0.0000 | -0.0000 0.0000 -0.0000 4.0000 0.0000 4.0000 -0.0000 0.0000 | ||
| 0 -0.0000 0.0000 -0.4000 0 0.4000 -0.0000 0.0000 | 0 -0.0000 0.0000 -0.4000 0 0.4000 -0.0000 0.0000 | ||
| + | |||
| + | |||
| >> k = 4; | >> k = 4; | ||
| >> u = cos(k*x) + 0.1*sin(k*x); | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| Line 115: | Line 104: | ||
| </code> | </code> | ||
| + | So the ordering of Fourier modes in the FFT for N=8 is k = [0 1 2 3 -4 -3 -2 -1] | ||
| + | --not what I suggested it would be in class. The normalization is also different: | ||
| + | the kth and -kth FFT coefficients for u(x) = cos(kx) are not 1/2, but N/2. | ||
| + | Also observe that the treatment of the highest-order Fourier mode k=-4 is a little weird, | ||
| + | in that the 8 represents the sum of the k=4 and k=-4 cosine modes, with no representation | ||
| + | of the k=4 and k=-4 sin modes. Without going into too much detail on this, Fourier spectral | ||
| + | PDE codes typically just zero the highest-order mode at every time step. We can ignore this | ||
| + | subtlety. | ||
| + | |||
| + | Based on the above, here's a short example of how to do Fourier differentiation in Matlab on | ||
| + | a [0,L] periodic domain. This should be plenty to get you started for the time-stepping codes. | ||
| + | |||
| + | <code> | ||
| + | % A demonstration of Fourier differentiation in Matlab | ||
| + | |||
| + | % Parameters | ||
| + | L = 10; | ||
| + | N = 16; | ||
| + | x = L/N*(0:N-1); | ||
| + | |||
| + | % Set up u(x) and du/dx for some known L-periodic function | ||
| + | a = 2*pi/L | ||
| + | u = sin(a*x) + 0.2*cos(a*2*x); | ||
| + | dudx = a*cos(a*x) - 0.4*a*sin(a*2*x); | ||
| + | |||
| + | % Now compute dudx via FFT | ||
| + | k = [0:N/2-1 -N/2:-1]; | ||
| + | alpha = 2*pi*k/L; | ||
| + | D = i*alpha; | ||
| + | dudx_fft = real(ifft(D.*fft(u))); | ||
| + | |||
| + | % Plot everybody | ||
| + | plot(x,u, 'b', x,dudx, 'g', x, dudx_fft, 'r.') | ||
| + | legend('u','du/dx','du/dx via FFT','location','southwest') | ||
| + | xlabel('x') | ||
| + | </code> | ||
gibson/teaching/spring-2014/iam950/hw1.1393515127.txt.gz · Last modified: 2014/02/27 07:32 by gibson
Page Tools
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 3.0 Unported
