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/25 07:52] gibson |
gibson:teaching:spring-2014:iam950:hw1 [2014/02/27 09:39] (current) gibson |
||
|---|---|---|---|
| Line 1: | Line 1: | ||
| ====== IAM 950 HW1 ====== | ====== IAM 950 HW1 ====== | ||
| - | **1.** Write numerical simulation for the heat eqn $u_t = \nu u_{xx}$ on a periodic domain | + | **Problem 1.** Write numerical simulation for the heat eqn $u_t = \nu u_{xx}$ on a periodic domain |
| - | $[-L_x/2, L_x/2]$ using Fourier spatial discretization and Crank-Nicolson temporal discretization. Verify that your code works correctly by simulating the Gaussian-decay solution | + | $[-L_x/2, L_x/2]$ using Fourier spatial discretization and Crank-Nicolson temporal discretization. |
| + | Verify that your code works correctly by simulating the Gaussian-decay solution | ||
| <latex> | <latex> | ||
| Line 11: | Line 12: | ||
| Initialize your numerical code with u(x,1) and then compare the results of the | Initialize your numerical code with u(x,1) and then compare the results of the | ||
| numerical time integration with u(x,t) evaluated from the above formula. I suggest plotting | numerical time integration with u(x,t) evaluated from the above formula. I suggest plotting | ||
| - | both quantities versus x at regular intervals in time. | + | both quantities versus x at regular intervals in time during the simulation, using "pause" (matlab) |
| + | to stop the simulation temporarily. | ||
| Make a plot of the numerical solution and the Gaussian solution versus x at t=100 to turn in. | Make a plot of the numerical solution and the Gaussian solution versus x at t=100 to turn in. | ||
| Are the two entirely consistent? If not, why not? | Are the two entirely consistent? If not, why not? | ||
| - | **2.** Copy and revise your heat equation code so that it simulated the 1d Swift-Hohenberg equation, | + | **Problem 2.** Copy and revise your heat equation code so that it simulated the 1d Swift-Hohenberg equation, |
| $u_t = (r-1) u - 2 u_{xx} - u_{xxxx} - u^3$ on a periodic domain $[0,L_x]$ using Fourier spatial | $u_t = (r-1) u - 2 u_{xx} - u_{xxxx} - u^3$ on a periodic domain $[0,L_x]$ using Fourier spatial | ||
| discretization and Crank-Nicolson/Adams-Bashforth semi-implicit temporal discretization. Use the | discretization and Crank-Nicolson/Adams-Bashforth semi-implicit temporal discretization. Use the | ||
| Line 22: | Line 24: | ||
| r = 0.2, Lx = 100, Nx = 256 (number of x gridpoints), dt = 1/16, and integrate from t=0 to 100. | r = 0.2, Lx = 100, Nx = 256 (number of x gridpoints), dt = 1/16, and integrate from t=0 to 100. | ||
| - | Make a colorplot of $u(x,t)$ in the $x,t$ plane (using every x gridpoint but plotting t at a larger interval than dt, perhaps at intervals of 1/2 or 1. It should look something like this | + | Make a colorplot of $u(x,t)$ in the $x,t$ plane (using every x gridpoint but plotting t at a larger |
| + | interval than dt, perhaps at intervals of 1/2 or 1. It should look something like this | ||
| {{:gibson:teaching:spring-2014:iam950:sh.png?direct&400|Simulations of Swift-Hohenberg eqn}} | {{:gibson:teaching:spring-2014:iam950:sh.png?direct&400|Simulations of Swift-Hohenberg eqn}} | ||
| + | |||
| + | When you are satisfied with the correctness of your simulation, make a plot of $u(x,100)$ versus x for | ||
| + | the above parameters and initial conditions to turn in. | ||
| Now play around with the forcing parameter in the range 0 < r < 1 and the initial condition, | Now play around with the forcing parameter in the range 0 < r < 1 and the initial condition, | ||
| Line 35: | Line 41: | ||
| - | **3.** Copy and revise your Swift-Hohenberg code so that it simulates the 1d Kuramoto-Sivashinsky equation, $u_t = - u_{xx} - u_{xxxx}- u u_x$, with parameters Lx = 100, Nx = 512, and dt = 1/16, | + | **Problem 3.** Copy and revise your Swift-Hohenberg code so that it simulates the 1d Kuramoto-Sivashinsky equation, $u_t = - u_{xx} - u_{xxxx}- u u_x$, with parameters Lx = 100, Nx = 512, and dt = 1/16, |
| initial condition $u(x,0) = \cos(ax) + 0.1 \cos(2ax)$ (same as before), and integrate from t=0 to t=100. Make a colorplot of $u(x,t)$ as before. It should look like this | initial condition $u(x,0) = \cos(ax) + 0.1 \cos(2ax)$ (same as before), and integrate from t=0 to t=100. Make a colorplot of $u(x,t)$ as before. It should look like this | ||
| {{:gibson:teaching:spring-2014:iam950:ks.png?direct&400|Kuramoto-Sivashinsky simulation}} | {{:gibson:teaching:spring-2014:iam950:ks.png?direct&400|Kuramoto-Sivashinsky simulation}} | ||
| - | Now play around with the domain size Lx (adjusting Nx so that Lx/Nx is approximately constant) and | + | Make a plot of $u(x,100)$ versus x to turn in. Now play around with the domain size Lx (adjusting |
| - | the initial condition, and answer the following questions | + | Nx so that Lx/Nx is approximately constant), try replacing the initial condition with random numbers |
| + | of order 0.01, and answer the following questions | ||
| - Does the wavelength of the final pattern depend on r or the initial condition? | - Does the wavelength of the final pattern depend on r or the initial condition? | ||
| Line 47: | Line 54: | ||
| - What qualitative changes in behavior do you observe as L increases from small (L=1) to large (L=100)? | - What qualitative changes in behavior do you observe as L increases from small (L=1) to large (L=100)? | ||
| - | It's probably a good idea to use small, completely random initial conditions as you vary L (so that all modes are excited initially in a way that's independent of L). | + | **Problem 4.** For each of the above codes, fix the initial condition and parameters, and run the simulation |
| + | 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. 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> | ||
| + | >> 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; | ||
| + | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| + | >> [real(fft(u)); imag(fft(u))] | ||
| + | ans = | ||
| + | 8 0 0 0 0 0 0 0 | ||
| + | 0 0 0 0 0 0 0 0 | ||
| + | |||
| + | |||
| + | >> k = 1; | ||
| + | >> u = cos(k*x) + 0.1*sin(k*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 | ||
| + | |||
| + | |||
| + | >> k = 2; | ||
| + | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| + | >> [real(fft(u)); imag(fft(u))] | ||
| + | ans = | ||
| + | -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 | ||
| + | |||
| + | |||
| + | >> k = 3; | ||
| + | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| + | >> [real(fft(u)); imag(fft(u))] | ||
| + | ans = | ||
| + | -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 | ||
| + | |||
| + | |||
| + | >> k = 4; | ||
| + | >> u = cos(k*x) + 0.1*sin(k*x); | ||
| + | >> [real(fft(u)); imag(fft(u))] | ||
| + | ans = | ||
| + | 0 0 0 0 8.0000 0 0 0 | ||
| + | 0 0.0000 0 0.0000 0 -0.0000 0 -0.0000 | ||
| + | </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.1393343578.txt.gz · Last modified: 2014/02/25 07:52 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
