AMATH 481/581 - Autumn 2022 Homework #5
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
AMATH 481/581 - Autumn 2022
Homework #5
Submissions accepted until 11:59 PM (PT) Wednesday, December 14, 2022
To submit this assignment, upload your main homework file (.m, .py, or .ipynb) to Grade- scope. Additionally you may upload a .pdf you create to demonstrate mastery of one or both presentation skills.
In this homework we will study a reaction-diffusion system. Reaction-diffusion equa- tions are used for modeling a variety of physical phenomena, especially in chemistry. Some examples include crystal growth and chemical reactions. I highly recommend look- ing up some of these examples, such as on Wikipedia or by looking up the “BZ reaction.” Reaction-diffusion equations are also very interesting mathematically because they can undergo a number of bifurcations leading to chaos (see AMATH 402 or 502 for more!).
A reaction-diffusion system for the two quantities U (z．y．t) and V (z．y．t), can be modeled by
Vt = A (a)U - λ(a)V + D2 V2 V． (2)
where a2 = U2 + V2 and V2 = θ北(2) + θy(2) . We will consider a particular system with
A (a) = -βa2． (4)
although you can try other choices of λ and A in your investigation of the system.
We will study this reaction-diffusion system with the spiral initial condition
U (z．y．0) = ╱tanh( ^z2 + y2 - a)、cos ╱a(z．y)m - ^z2 + y2 、 (5)
V (z．y．0) = ╱tanh( ^z2 + y2 - a)、sin ╱a(z．y)m - ^z2 + y2、． (6)
where a(z．y) is the angle between the z axis and the complex number z + iy, and a and m are numbers to be given later. The function a is implemented in MATLAB using angle(x+1i*y) and in python using np .angle(x+1j*y). See the Canvas page for code
implementing these initial conditions in MATLAB and python.
We will solve this initial-value problem using a Fourier spectral method and a Cheby- shev pseudo-spectral method, for different boundary conditions. For periodic boundary conditions you will use the spectral FFT method and for Dirichlet boundary conditions you will use the Chebyshev method. Note that we have two coupled ODEs in two di- mensions, so we need to use a reshaped meshgrid and also stack U and V to create one solution vector. In both cases we will solve using RK45.
In both of our solves we will use z．y e [-10．10], β = 1, D1 = D2 = 0．1, and will solve for t e [0．25] with a step size of ∆t = 0．5. In the end I highly recommend animating your solutions using pcolor, because the result is very cool.
1. We will first solve the equation on the periodic domain using the FFT spectral method. To do so, set m = 3, a = 0, and use n = 64 points in the z and y directions (making sure not to include endpoints twice). When creating the k vector, we now do not need to (and should not) set the first value of k to 10 −6 because we do not need to divide by k .
Deliverables:
(a) Create the mesh for the z and y variables. Save the Ⅹ mesh to the variable A1 to verify that you have setup the grid correctly. X should be a 64 × 64 matrix.
(b) Define the initial condition for U (5) on the grid and save it to the variable
A2. A2 should be a 64 × 64 matrix.
(c) Take the FFT of the initial conditions using fft2, call them 
0 and 
0 . Note that these will have a real and imaginary part. Save the real part of 
0 (real MATLAB and np .real python) of 
0 to the variable A3. A3 should be a
64 × 64 matrix.
(d) In order to solve using RK45 we need the initial condition to be a vector/array. Reshape 
0 and 
0 into arrays for U and V and then stack them into a column vector of dimension 8192 × 1 (note, 8192 = 64 × 64 × 2). This whole vector will also have a real and imaginary part. Save the imaginary part (imag MATLAB and np .imag python) of the initial condition to the variable A4.
(e) Timestep using RK45 to get the solution for t e [0．25] and save the solution with ∆t = 0．5. Note that this is Fourier Transform of the solution at the 51 time points. Save the real part of the result as a 8192 × 51 matrix to the variable A5. Save the imaginary part of the result as a 8192 × 51 matrix to the variable A6.
(f) Save the real part of the Fourier transform of the solution U at t = 2 as a 4096 × 1 vector to the variable A7.
(g) Reshape the Fourier transform of U at t = 2 to a 64 × 64 matrix. Save the real part to the variable A8. Make sure that the way you are reshaping satisfies, e.g., if
: = ┌ 1 2 3 4 5 6 7 8 9]．
then the reshaped 3 × 3 matrix should be
┌ 1 4
'
'
'3 6
7┐
'
'
9'
(h) Save the inverse Fourier transform of the above result (using ifft2, giving U (z．y．2)), a 64 × 64 matrix, to the variable A9. This is the quantity you want to plot to see cool stuff .
2. You will now solve the same problem with Dirichlet Boundary conditions,
U (z = -10．y．t) = U (z = 10．y．t) = U (z．y = -10．t) = U (z．y = 10．t) = 0．
using the Chebyshev pseudo-spectral method. To do so, set m = 2, a = 1, and use n = 30 when creating the Chebyshev matrix and Chebyshev points. You should use the functions I provided for creating the matrix and points.
Deliverables
(a) Create the Chebyshev differentiation matrix for the Laplacian V2 = θx(2) + θy(2) . Scale it accordingly to match our domain and remove the first and last rows and columns to account for the boundary conditions. Save the resulting 841 × 841 matrix to the variable A10.
(b) Create the Chebyshev points in the z and y directions with the correct scaling for our domain and with the boundary points removed. Create the mesh and
save the Y mesh, a 29 × 29 matrix, to the variable A11.
(c) Define the initial condition for V (6) on the grid and save it to the variable A12. This will be a 29 × 29 matrix.
(d) In order to solve using RK45 we need the initial condition to be a vector/array. Reshape the initial conditions for U and V and stack them (U on top of V) to create a 1682 × 1 vector. Save it to the variable A13. (Note, 1682 = 29×29×2).
(e) Timestep using RK45 to get the solution for t e [0．25] and save the solution
with ∆t = 0．5. Save the result as a 51 × 1682 matrix to the variable A14. (f) Save V at time t = 2 as a 841 × 1 vector to the variable A15.
(g) Reshape U(z．y．t = 2) and V (z．y．t = 2) (at t = 2) to 29 × 29 matrices. Make sure to reshape these according the specifications in the problem above. After reshaping, add in the boundary conditions at z = -10 and z = 10 by padding the matrix with 0s in the first and last rows and columns. Save the resulting 31 × 31 matrix for V at t = 2 to the variable A16. You may want to animate this quantity V (z．y．t).
Presentation skills: Use pcolor to animate U(z．y．t) on the periodic domain. Be sure to remove the grid, label your axes, include a color bar, use a constant axis limits for the colors (to clearly see the behavior), and a title name that changes in time, showing at which time the solution is being shown. Your animation should not be choppy. If it looks too choppy, decrease ∆t and try again. Submit your animation to Canvas under
Homework 5.
Other fun things, optional: Besides these answers, you may want to investigate and construct various one- and two-armed (m = 1．2) spirals for this system. Also, investigate when the solutions become unstable and “chaotic” in nature. Investigate the system for all three (Dirichlet, Neumann, and periodic) boundary conditions. You may also want to play around with increasing the amount of diffusion (D1 and D2 ) as well as the amount of advection (β).
2022-12-14