-# Introduction\r
-\r
-This program is used to solve the generalized Swift-Hohenberg equation\r
-\r
-$$\begin{aligned}\r
- \frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3\r
-\end{aligned}$$\r
-\r
-where $k_c$ is the wave number, $r$ is some fixed constant, and\r
-$g_1$ is a parameter which determines the behavior of the solutions.\r
-Note that the equation is simply called the Swift-Hoheneberg equation if\r
-$g_1 = 0$. For this solver, we chose $k_c = 1$ and $r = 0.3$.\r
-Choosing $k_c$ to be 1 will mean that our solutions have a pattern\r
-wavelength of $2\pi$. We choose $r = 0.3$ because solutions are\r
-reasonably well behaved for small values of $r$ and $g_1$, but there\r
-are interesting behaviors that occur when $g_1$ is smaller or larger\r
-than $r$ in magnitude, so this allows us room to vary $g_1$ and\r
-explore these behavior. To summarize, this code solves:\r
-\r
-$$\begin{aligned}\r
- \frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3\r
-\end{aligned}$$\r
-\r
-# Discretization and Solving the Bilaplacian\r
-\r
-The equation has two aspects which are challenging to solve: the\r
-nonlinear terms $g_1u^2 - u^3$ and the Bilaplacian operator\r
-$(1 + \Delta)^2$, which introduces $4^{th}$ derivatives. To deal\r
-with the Bilaplacian, we introduce a variable $v$ and construct a\r
-system of PDEs:\r
-\r
-$$\begin{aligned}\r
- \frac{\partial u}{\partial t} &= 0.3u - (1 + \Delta) v + g_1 u^2 - u^3\\\r
- (1 + \Delta)u &= v\r
-\end{aligned}$$\r
-\r
-We can solve these two equations simultaneously by treating our\r
-finite elements as vector valued, and interpreting our system of\r
-equations as a single vector-valued PDE. We can handle the nonlinear\r
-terms by treating them fully explicitly. If we discretize in time and\r
-rearrange terms, our system of equations becomes \r
-\r
-$$\begin{aligned}\r
- (1 - kr)U_n + k(1 + \Delta)V_n &= U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
- (1 + \Delta)U_n - V_n &= 0\r
-\end{aligned}$$ \r
-\r
-where $k$ is the discrete timestep, $U_n$ and\r
-$V_n$ are the solutions for $u$ and $v$ at the current timestep,\r
-and $U_{n-1}$ and $V_{n-1}$ are the solutions for $u$ and $v$ at\r
-the previous timestep. We then reframe this system as a vector valued\r
-problem\r
-\r
-$$\begin{aligned}\r
- \left(\begin{matrix}\r
- 1 - kr & k(1 + \Delta)\\\r
- 1 + \Delta & -1\r
- \end{matrix}\right)\r
- \left(\begin{matrix}\r
- U_n\\\r
- V_n\r
- \end{matrix}\right) &= \left(\begin{matrix}\r
- U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
- 0\r
- \end{matrix}\right)\r
-\end{aligned}$$\r
-\r
-As usual, we multiply each side of the equation by a\r
-test function \r
-\r
-$$\overrightarrow\varphi_i = \left(\begin{matrix}\r
- \phi_i\\\r
- \psi_i\r
-\end{matrix}\right)$$\r
-\r
-and then integrate over the domain $\Omega$ to get the equation\r
-\r
-$$\begin{aligned}\r
- \int_\Omega \left(\begin{matrix}\r
- \phi_i\\ \r
- \psi_i\r
- \end{matrix}\right)\cdot\left(\begin{matrix}\r
- 1 - kr & k(1 + \Delta)\\\r
- 1 + \Delta & -1\r
- \end{matrix}\right)\r
- \left(\begin{matrix}\r
- U_n\\\r
- V_n\r
- \end{matrix}\right) &= \int_\Omega \left(\begin{matrix}\r
- \phi_i\\\r
- \psi_i\r
- \end{matrix}\right)\cdot\left(\begin{matrix}\r
- U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
- 0\r
- \end{matrix}\right)\\\r
-\end{aligned}$$\r
-\r
-We can expand our solution vector in this basis\r
-\r
-$$\begin{aligned}\r
- \int_\Omega \sum_j u_j\left(\begin{matrix}\r
- \phi_i\\\r
- \psi_i\r
- \end{matrix}\right)\cdot\left(\begin{matrix}\r
- 1 - kr & k(1 + \Delta)\\\r
- 1 + \Delta & -1\r
- \end{matrix}\right)\r
- \left(\begin{matrix}\r
- \phi_j\\\r
- \psi_j\r
- \end{matrix}\right) &= \int_\Omega\left(\begin{matrix}\r
- \phi_i\\\r
- \psi_i\r
- \end{matrix}\right)\cdot\left(\begin{matrix}\r
- U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\\r
- 0\r
- \end{matrix}\right)\r
-\end{aligned}$$\r
-\r
-and finally expand out the matrix multiplication\r
-and dot products, then apply the divergence theorem to obtain\r
-a single equation:\r
-\r
-$$\begin{aligned}\r
- \sum_j u_j \int_\Omega[(1 - kr)\phi_i\phi_j + k\phi_i\psi_j - k\nabla\phi_i\nabla\psi_j + \psi_i\phi_j - \nabla\psi_i\nabla\psi_j - \psi_i\psi_j] &= \int_\Omega\phi_i(U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3)\r
-\end{aligned}$$\r
-\r
-This last equation represents matrix multiplication of the\r
-solution vector by the $i^{th}$ row of the system matrix, and the left\r
-hand side without the summation or $u_j$ term is what we use to\r
-compute the $(i, j)^{th}$ entry of the system matrix.\r
-\r
-# Boundary Conditions and Choosing a Suitable Domain\r
-\r
-This code implements both zero Dirichlet and zero Neumann boundary\r
-conditions. Both of these conditions affect the patterns that form. To\r
-minimize this effect, we use Neumann boundary conditions and we choose\r
-the boundary to be some multiple of the wave number. For example, this\r
-code chooses the square mesh to have a side length of $6\cdot 2\pi$.\r
-For all other domains used, we chose a domain size with a similar area\r
-to that of the square. For instance, we solve on a torus with outer\r
-radius 9 and inner radius 4 because this results in exactly the same\r
-area as the square. Note that this is not strictly necessary for the\r
-code to function, but does make it easier to compare results between\r
-different geometries.\r
-\r
-# Initial Conditions and Parameters\r
-\r
-The code implements two main types of initial conditions: random initial\r
-conditions, and creating a small initial hot spot. The SH equation is\r
-interesting because it describes pattern formation and\r
-self-organization, so choosing random initial conditions allows for this\r
-to be observed. Note that the results shown below were all run with the\r
-initial seed 314, which was arbitrarily chosen. Setting a fixed seed is\r
-useful for comparing pattern formation with different choices of\r
-parameters in the SH equation.\r
-\r
-The hot spot initial condition is useful for the opposite reason: it is\r
-very simple, but it lets us see what happens to a single pattern\r
-\"wave\" as it propagates along our surface. This is particularly\r
-useful in distinguishing the effect of curvature and geometry on pattern\r
-propagation.\r
-\r
-As previously mentioned, we chose $k_c = 1$ and $r = 0.3$ for this\r
-code. We then let $g_1$ be the parameter that we change to vary the\r
-patterns formed. On the plane, increasing the value of $g_1$ allows\r
-for the formation of hexagonal grids rather than just ripples. Varying\r
-$g_1$ does something similar to patterns on a curved manifold, though\r
-with notably different effects in some cases. Increasing $g_1$ also\r
-causes the solution to grow larger in magnitude at certain points.\r
-\r
-# Checking Convergence\r
-\r
-We checked the convergence of this code using 3 tests: we confirmed that\r
-a constant initial condition remained constant and converged to a\r
-solution that was verified using an ordinary differential equation, we\r
-checked that solutions on the square converged across mesh refinements,\r
-and we checked that solutions converged over refinements of the timestep\r
-on the finest mesh.\r
-\r
-Below are the results of several runs of constant initial conditions\r
-\r
-\r
-\r
-We also validated that given a fixed random start on a very fine mesh,\r
-refining the timestep resulted in the same final solution. The initial\r
-condition for each is shown above, While the final solutions are shown in the matrix below. Note that the\r
-timestep begins at 1/25 and the denominator increases by 25 across each\r
-row, to a max of 1/200 in the bottom right:\r
-\r
-\r
-\r
-We validated that solutions converged across mesh refinement by defining\r
-psuedorandom functions\r
-$\displaystyle f(x) = \sum_{n=1}^{10} C_n \sin\left(\frac{16x}{3i}\right)$\r
-and\r
-$\displaystyle g(y) = \sum_{n=1}^{10} D_n \sin\left(\frac{16y}{3i}\right)$,\r
-where $C_i$ and $D_i$ are randomly chosen in the range\r
-$(-\sqrt{r}, \sqrt{r})$. The overall pseudorandom function is\r
-$h(x) = f(x)g(y)$. Note that the period of the sine waves was chosen\r
-so that the smallest wave could be resolved by a mesh refinement of 7 or\r
-higher. The following matrix shows the initial and final solution\r
-ranging from a refinement of 0 to a refinement of 7:\r
-\r
-\r
-\r
-\r
-\r
-# Results\r
-\r
-We can see the effects of varying the $g_1$ parameter and the effects\r
-of curvature using the hot spot initial condition. On the plane, an\r
-initial hot spot creates one ripple wave, which breaks into discrete\r
-pieces as $g_1$ is increased. In the matrix below, $g_1$ is\r
-increased by 0.2 starting from 0 to a maximum value of 1.4. Note that\r
-each final solution is at 100 time units:\r
-\r
-\r
-\r
-On the cylinder, the front looks similar to the square, but the back has\r
-an overlapping wave pattern:\r
-\r
-\r
-\r
-\r
-On the sphere, the hot spot generates a single wave. Note that this may\r
-be due to the fact that our sphere has a surface area proportional to\r
-the period of our pattern wave.\r
-\r
-\r
-\r
-On the torus, the pattern propagates similar to the cylinder, with some\r
-minor imperfections\r
-\r
-\r
-\r
-But on the back side of the torus, we see wave overlapping and spot\r
-patterns forming\r
-\r
-\r
-\r
-On shapes with stranger curvature, we can see that the pattern wave has\r
-a tendency to break apart when crossing lines of curvature. This shape\r
-was made by warping the boundary of a cylinder by a cosine wave, and is\r
-equivalent to the surface of revolution bounded by\r
-$1 + 0.5\cos(\frac{\pi}{10}x)$\r
-\r
-\r
-\r
-\r
-\r
-Finally, here is a small selection of random initial conditions and the\r
-patterns that form. Each image sequence was taken at times 0, 10, 25,\r
-50, and 100:\r
-\r
-\r
-\r
-\r
-\r
+# Introduction
+
+This program is used to solve the generalized Swift-Hohenberg equation
+
+$$\begin{aligned}
+ \frac{\partial u}{\partial t} = ru - (k_c + \Delta)^2 u + g_1 u^2 - u^3
+\end{aligned}$$
+
+where $k_c$ is the wave number, $r$ is some fixed constant, and
+$g_1$ is a parameter which determines the behavior of the solutions.
+Note that the equation is simply called the Swift-Hoheneberg equation if
+$g_1 = 0$. For this solver, we chose $k_c = 1$ and $r = 0.3$.
+Choosing $k_c$ to be 1 will mean that our solutions have a pattern
+wavelength of $2\pi$. We choose $r = 0.3$ because solutions are
+reasonably well behaved for small values of $r$ and $g_1$, but there
+are interesting behaviors that occur when $g_1$ is smaller or larger
+than $r$ in magnitude, so this allows us room to vary $g_1$ and
+explore these behavior. To summarize, this code solves:
+
+$$\begin{aligned}
+ \frac{\partial u}{\partial t} = 0.3u - (1 + \Delta)^2 u + g_1 u^2 - u^3
+\end{aligned}$$
+
+# Discretization and Solving the Bilaplacian
+
+The equation has two aspects which are challenging to solve: the
+nonlinear terms $g_1u^2 - u^3$ and the Bilaplacian operator
+$(1 + \Delta)^2$, which introduces $4^{th}$ derivatives. To deal
+with the Bilaplacian, we introduce a variable $v$ and construct a
+system of PDEs:
+
+$$\begin{aligned}
+ \frac{\partial u}{\partial t} &= 0.3u - (1 + \Delta) v + g_1 u^2 - u^3\\
+ (1 + \Delta)u &= v
+\end{aligned}$$
+
+We can solve these two equations simultaneously by treating our
+finite elements as vector valued, and interpreting our system of
+equations as a single vector-valued PDE. We can handle the nonlinear
+terms by treating them fully explicitly. If we discretize in time and
+rearrange terms, our system of equations becomes
+
+$$\begin{aligned}
+ (1 - kr)U_n + k(1 + \Delta)V_n &= U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\
+ (1 + \Delta)U_n - V_n &= 0
+\end{aligned}$$
+
+where $k$ is the discrete timestep, $U_n$ and
+$V_n$ are the solutions for $u$ and $v$ at the current timestep,
+and $U_{n-1}$ and $V_{n-1}$ are the solutions for $u$ and $v$ at
+the previous timestep. We then reframe this system as a vector valued
+problem
+
+$$\begin{aligned}
+ \left(\begin{matrix}
+ 1 - kr & k(1 + \Delta)\\
+ 1 + \Delta & -1
+ \end{matrix}\right)
+ \left(\begin{matrix}
+ U_n\\
+ V_n
+ \end{matrix}\right) &= \left(\begin{matrix}
+ U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\
+ 0
+ \end{matrix}\right)
+\end{aligned}$$
+
+As usual, we multiply each side of the equation by a
+test function
+
+$$\overrightarrow\varphi_i = \left(\begin{matrix}
+ \phi_i\\
+ \psi_i
+\end{matrix}\right)$$
+
+and then integrate over the domain $\Omega$ to get the equation
+
+$$\begin{aligned}
+ \int_\Omega \left(\begin{matrix}
+ \phi_i\\
+ \psi_i
+ \end{matrix}\right)\cdot\left(\begin{matrix}
+ 1 - kr & k(1 + \Delta)\\
+ 1 + \Delta & -1
+ \end{matrix}\right)
+ \left(\begin{matrix}
+ U_n\\
+ V_n
+ \end{matrix}\right) &= \int_\Omega \left(\begin{matrix}
+ \phi_i\\
+ \psi_i
+ \end{matrix}\right)\cdot\left(\begin{matrix}
+ U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\
+ 0
+ \end{matrix}\right)\\
+\end{aligned}$$
+
+We can expand our solution vector in this basis
+
+$$\begin{aligned}
+ \int_\Omega \sum_j u_j\left(\begin{matrix}
+ \phi_i\\
+ \psi_i
+ \end{matrix}\right)\cdot\left(\begin{matrix}
+ 1 - kr & k(1 + \Delta)\\
+ 1 + \Delta & -1
+ \end{matrix}\right)
+ \left(\begin{matrix}
+ \phi_j\\
+ \psi_j
+ \end{matrix}\right) &= \int_\Omega\left(\begin{matrix}
+ \phi_i\\
+ \psi_i
+ \end{matrix}\right)\cdot\left(\begin{matrix}
+ U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3\\
+ 0
+ \end{matrix}\right)
+\end{aligned}$$
+
+and finally expand out the matrix multiplication
+and dot products, then apply the divergence theorem to obtain
+a single equation:
+
+$$\begin{aligned}
+ \sum_j u_j \int_\Omega[(1 - kr)\phi_i\phi_j + k\phi_i\psi_j - k\nabla\phi_i\nabla\psi_j + \psi_i\phi_j - \nabla\psi_i\nabla\psi_j - \psi_i\psi_j] &= \int_\Omega\phi_i(U_{n-1} + kg_1U_{n-1}^2 - kU_{n-1}^3)
+\end{aligned}$$
+
+This last equation represents matrix multiplication of the
+solution vector by the $i^{th}$ row of the system matrix, and the left
+hand side without the summation or $u_j$ term is what we use to
+compute the $(i, j)^{th}$ entry of the system matrix.
+
+# Boundary Conditions and Choosing a Suitable Domain
+
+This code implements both zero Dirichlet and zero Neumann boundary
+conditions. Both of these conditions affect the patterns that form. To
+minimize this effect, we use Neumann boundary conditions and we choose
+the boundary to be some multiple of the wave number. For example, this
+code chooses the square mesh to have a side length of $6\cdot 2\pi$.
+For all other domains used, we chose a domain size with a similar area
+to that of the square. For instance, we solve on a torus with outer
+radius 9 and inner radius 4 because this results in exactly the same
+area as the square. Note that this is not strictly necessary for the
+code to function, but does make it easier to compare results between
+different geometries.
+
+# Initial Conditions and Parameters
+
+The code implements two main types of initial conditions: random initial
+conditions, and creating a small initial hot spot. The SH equation is
+interesting because it describes pattern formation and
+self-organization, so choosing random initial conditions allows for this
+to be observed. Note that the results shown below were all run with the
+initial seed 314, which was arbitrarily chosen. Setting a fixed seed is
+useful for comparing pattern formation with different choices of
+parameters in the SH equation.
+
+The hot spot initial condition is useful for the opposite reason: it is
+very simple, but it lets us see what happens to a single pattern
+\"wave\" as it propagates along our surface. This is particularly
+useful in distinguishing the effect of curvature and geometry on pattern
+propagation.
+
+As previously mentioned, we chose $k_c = 1$ and $r = 0.3$ for this
+code. We then let $g_1$ be the parameter that we change to vary the
+patterns formed. On the plane, increasing the value of $g_1$ allows
+for the formation of hexagonal grids rather than just ripples. Varying
+$g_1$ does something similar to patterns on a curved manifold, though
+with notably different effects in some cases. Increasing $g_1$ also
+causes the solution to grow larger in magnitude at certain points.
+
+# Checking Convergence
+
+We checked the convergence of this code using 3 tests: we confirmed that
+a constant initial condition remained constant and converged to a
+solution that was verified using an ordinary differential equation, we
+checked that solutions on the square converged across mesh refinements,
+and we checked that solutions converged over refinements of the timestep
+on the finest mesh.
+
+Below are the results of several runs of constant initial conditions
+
+
+
+We also validated that given a fixed random start on a very fine mesh,
+refining the timestep resulted in the same final solution. The initial
+condition for each is shown above, While the final solutions are shown in the matrix below. Note that the
+timestep begins at 1/25 and the denominator increases by 25 across each
+row, to a max of 1/200 in the bottom right:
+
+
+
+We validated that solutions converged across mesh refinement by defining
+psuedorandom functions
+$\displaystyle f(x) = \sum_{n=1}^{10} C_n \sin\left(\frac{16x}{3i}\right)$
+and
+$\displaystyle g(y) = \sum_{n=1}^{10} D_n \sin\left(\frac{16y}{3i}\right)$,
+where $C_i$ and $D_i$ are randomly chosen in the range
+$(-\sqrt{r}, \sqrt{r})$. The overall pseudorandom function is
+$h(x) = f(x)g(y)$. Note that the period of the sine waves was chosen
+so that the smallest wave could be resolved by a mesh refinement of 7 or
+higher. The following matrix shows the initial and final solution
+ranging from a refinement of 0 to a refinement of 7:
+
+
+
+
+
+# Results
+
+We can see the effects of varying the $g_1$ parameter and the effects
+of curvature using the hot spot initial condition. On the plane, an
+initial hot spot creates one ripple wave, which breaks into discrete
+pieces as $g_1$ is increased. In the matrix below, $g_1$ is
+increased by 0.2 starting from 0 to a maximum value of 1.4. Note that
+each final solution is at 100 time units:
+
+
+
+On the cylinder, the front looks similar to the square, but the back has
+an overlapping wave pattern:
+
+
+
+
+On the sphere, the hot spot generates a single wave. Note that this may
+be due to the fact that our sphere has a surface area proportional to
+the period of our pattern wave.
+
+
+
+On the torus, the pattern propagates similar to the cylinder, with some
+minor imperfections
+
+
+
+But on the back side of the torus, we see wave overlapping and spot
+patterns forming
+
+
+
+On shapes with stranger curvature, we can see that the pattern wave has
+a tendency to break apart when crossing lines of curvature. This shape
+was made by warping the boundary of a cylinder by a cosine wave, and is
+equivalent to the surface of revolution bounded by
+$1 + 0.5\cos(\frac{\pi}{10}x)$
+
+
+
+
+
+Finally, here is a small selection of random initial conditions and the
+patterns that form. Each image sequence was taken at times 0, 10, 25,
+50, and 100:
+
+
+
+
+
\ No newline at end of file
projects / code-gallery.git / commitdiff