Mixed Finite Element Formulation and the necessity $H(\textrm{div})$

One can solve Poisson’s problem $-\Delta u = f$ in $d$ dimensions with homogeneous Dirichlet boundary conditions using a mixed formulation as explained below:

Let $\sigma = \nabla u$, then for a sufficiently smooth function $\tau$, by Green’s theorem
\begin{align*}
(\sigma, \tau) &= (\nabla u, \tau) \\
&= -(u, \textrm{div } \tau).
\end{align*}
Again, choosing $v$ a function sufficiently smooth, we have
\begin{align*}
f = -\textrm{div } \sigma \implies (f, v) = (-\textrm{div } \sigma, v).
\end{align*}
This gives the saddle-point problem: find $(u, \sigma) \in V \times M$ such that
\begin{align*}
(\sigma, \tau) + (u, \textrm{div } \tau) &= 0\\
(\textrm{div } \sigma, v) &= -(f, v)
\end{align*}
hold for all $(\tau, v) \in V \times M$. Note that we don’t have to take a derivative of $u$, hence it’s natural to try $M = L^2$, but what about the space $V$?

One very easy choice to guess is $V = [H^1(\Omega)]^d$ as we want the divergence to be all defined, but unfortunately this doesn’t work as the gradient of the solution to Poisson’s problem can easily not be in $[H^1(\Omega)]^d$.

In order to illustrate this, consider $u =\left(r^{2/3}-r^{5/3}\right) \sin \left(\frac{2 \theta }{3}\right)$ on the domain of the unit circle with bottom left quarter taken out. It’s not hard to see that $u = 0$ on the boundary of the domain, and we can easily find the $f$ such that it satisfies Poisson’s equation. Now, we can either calculate the gradient exactly or argue as follows.

First, recall how to take a gradient in polar coordinates. Note that $\partial_r u \approx r^{-1/3}$ plus higher order terms and that $\frac{1}{r}\partial_\theta u \approx r^{-1/3}$ plus higher order terms also. Now, one can easily calculate the $H^1$ seminorm to see that the derivative is unbounded as we’re integrating over $[0,1]$ with $(r^{-4/3})^2r$ terms (the extra $r$ comes from the change of basis from polar integration).

The above is an example of why the space $H(\textrm{div})$ is needed.

Celeste

Once again, Dunkey has proven himself to be a modern day Donald Draper… in some sense. I bought Celeste almost strictly due to how fun it seemed. It truly is a great game with tight controls and extremely interesting level design.

The first point is starting to get quite standard now, but I want to reiterate on the latter point. It seems many high-ceiling platformers like Supermeatboy or the age-old N game constantly rely on precision. Celeste throws that away with an emphasis on when/where/how you use the air dash. That air dash, that one extra mechanic, really is the crux. Honestly, it reminded me of Ori’s dash, but the level design here is more like a puzzle.

Of course, the game can get a bit annoying. There are places where precise timing is the only way through the level (or so it seems to me?). Flag number 9 is particularly annoying. One can turn on the assist mode, but it really breaks the game by making it very easy. Another quirk is that it doesn’t save automatically when quitting from the Switch; this caused me to have to beat certain levels twice as I was switching between games.

All in all, quite a fun game with a ton of content. Worth it for just $20.

A note on a paper by Babuska

This is a post regarding the paper Efficient Preconditioning for the $p$-Version Finite Element Method in Two Dimensions. In Lemma 3.3, the basis required is the bubble (or interior) polynomials, the set of edge functions orthogonal to the interiors, and linear (of bilinear) functions.

The edge functions orthogonality is explicitly stated as $\hat a(u, v) = 0$ for all $u \in \Gamma_i$ (edge space) and $v \in \mathcal{I}$ interior space. A very natural question is why the vertex functions does not need to be orthogonal to the interior functions. The fun fact is that it secretly is.

Note that the paper is for the $H^1$ semi-norm, hence $\hat a(u, v) = \int_T \nabla u \nabla v \, dx$. Now let $u$ be a hat (or bilinear) function, and let $v \in \mathcal{I}$. Then we have that
\begin{align*}
\int_T \nabla u \nabla v = \int_{\partial T} u \partial_n v – \int_T \Delta u v = 0.
\end{align*}
The first term is 0 due to the bubble functions vanishing on the boundary, and that $\Delta u = 0$ because it is linear.

Gray-Scott Equations

For our first two papers, we essentially reused the same few examples as model problems to test our method with (sine-Gordon and Brusselator). For our next paper, my advisor wanted something different and pointed towards the Grey-Scott equations. It’s a simple reaction diffusion equation as follows
\begin{align*}
\frac{\partial u}{\partial t} &= d_u\Delta u – uv^2 + F(1 – u) \\
\frac{\partial v}{\partial t} &= d_v\Delta v + uv^2 – (F+k)v
\end{align*}
where $F, k, d_u, d_v$ are constants.

There’s a short paper (“Complex Patterns in a Simple System,” by John E. Pearson) where he plots the function for different values of $F, k$. The problem for me in replicating that paper is that Pearson employed a periodic boundary condition, which is easy to implement for finite different and spectral methods, but a bit awkward for finite element methods (if you’re not using a very specific mesh).

The solution is quite nifty. Instead of a 2D plane, we simply project the domain onto a torus. It turns out FEM code on a surface is almost the same as on a plane, unless one uses curvilinear elements which then becomes a hassle. Furthermore, visualization gives pretty cool results… take a look at the two simulations below (their parameters differ just ever so slightly, but ends up giving drastically different patterns, though both reach a steady state):

New Year New Notes

Every year, I make the resolution of blogging and journalling more, and each passing year, I’ve noticed my failure. Yet here I am again stubborn as always to see how I can keep this up. Since it’s been awhile, here’s a laundry list of updates:

  1. Family and I went cruising on RCI to Haiti, Jamaica and Mexico over Christmas, and it felt just so artificial. There was a thin veneer of luxury painted over the backs of the employees, who were (almost) all from poor countries. Everything seems to be nickel and dimed…. nevertheless, it’s still relaxing not to do any chores or cooking.
  2. It’s incredibly hard to attempt to work at a place where you have never worked before. I look forward to getting back to writing and math tomorrow (and lifting/working out).
  3. Finished Sapiens, Fahrenheit 451 and How I Killed Pluto and Why It Had it Coming. All fun books to read, with Sapiens being a bit too long and Fahrenheit too short. The Pluto book was quite fun to read, and it’s interesting to see how other fields work.

P&B #155

Let $a, b, c$ be real numbers. Prove that three roots of the equation
\begin{align*}
\frac{b+c}{x-a} + \frac{c+a}{x-b} + \frac{a+b}{x-c} = 3
\end{align*}
are real.

Solution: The solution which the book lists is actually much more elegant than what I will present here. We let $s = a+b+c$ to be the sum. Inserting this into the equation, we see that
\begin{align*}
\frac{s-a}{x-a} + \frac{s-b}{x-b} + \frac{s-c}{x-c} = 3
\end{align*}
Now, multiplying both sides by $(x-a)(x-b)(x-c)$ and rearranging results in
\begin{align*}
(s-a)(x-b)(x-c) + (s-b)(x-a)(x-c) + \\
(s-c)(x-a)(x-b) – 3(x-a)(x-b)(x-c) = 0.
\end{align*}
By inspection, it’s pretty clear that $x = s$ is a root, hence if we find one additional real root we are done as imaginary roots must exist in pairs.

Here is where I deviated against the solutions. I plugged in $x = a, b, c$ respectively and considered all the cases separated by magnitude and showed that there must exist a pair to the left/right of $s$ such that it is positive and negative, hence there exists another real root.

The Pattern is a Lie

We started a new problem recently, but still within the realms of preconditioning. The exact question and algorithms will remain under wraps right now, but the preliminary results I got were incredibly curious; here’s a plot of condition number below:

Note that a constant condition number as we increase order is extremely good, but my adviser and I were not rejoicing.  It’s just that the result was too good! We were not expecting the condition number to not increase, as that’s far better than what current literature (and our intuition) suggests.

After a week of fumbling around, I’ve extended the computation and produced the following plot:

Not that the above plot is the same algorithm, just with extended order (or p). And ahhh, there it is,  the growth that we were looking for.

Even though it’s a worse result that what we initially saw, we can actually prove this version!

P&B #106

This is a problem from Putnam and Beyond:

Given that $9a^2 + 8ab + 7b^2 \le 6$, show that $7a + 5b + 12ab \le 9$.

My solution is extremely janky; nevertheless, it is still a solution that is slightly different from the book’s.

Our hypothesis implies that $-12ab \ge -9 + \frac{27}{2}a^2 + \frac{21}{2}b^2$ by manipulation. Note that our conclusion is equivalent to $9 – 7a – 5b – 12ab \ge 0$. Hence, we have

\begin{align*}
9 – 7a – 5b – 12ab &\ge \frac{27}{2}a^2 + \frac{21}{2}b^2- 7a – 5b \stackrel{?}{\ge} 0.
\end{align*}

Now consider the two quadratics $\frac{27}{2}a^2 – 7a$ and $\frac{21}{2}b^2 – 5b$. They are negative between $0$ and their other root, which is $a = \frac{14}{27}$ and $b = \frac{10}{21}$. Since the objective function is strictly increasing for in the positive numbers, checking the last root reveals that it is indeed satisfied as $7a + 5b + 12 ab = \frac{1696}{189} < 9$. Again, janky but doable in a reasonable time. (Of course, there's more proper methods with Lagrange multipliers and such...)

A Geometry Problem

Here’s an excellent question from AoPS:

In the figure to the left, circle $B$ is tangent to circle $A$ at $X$, circle $C$ is tangent to circle $A$ at $Y$, and circles $B$ and $C$ are tangent to each other. If $AB = 6, AC = 5,$ and $BC = 9$, what is $AX$?

The solution really requires no geometry besides the very simple definitions of a circle and a little bit of algebra. Let $r_B, r_C$ be the radius of circles $B, C$ respectively, then we have $r_B + r_C = BC = 9$. Second, $AX = AY = AC + CY = 5 + r_C = 6 + r_B$. Solving these two systems of equation gives 4 and 5. Hence the radius is 10.