Parabolic finite-element methods — Examples

Here $a(x,y)=xy$ has $M=\alpha =1$. The derivative is $u^{\mathrm{\prime }}(t)=\frac{3}{4}{t}^{-1\mathrm{/}4}$ and ${\int }_0^1\mathrm{\mid }{u}^{\mathrm{\prime }}(t){\mathrm{\mid }}^{2}dt=\frac{9}{16}{\int }_0^1{t}^{-1\mathrm{/}2}dt=\frac{9}{8}<\mathrm{\infty }$; also $u\in L^2(0,1)$, so $u\in H^1(0,1)$ and $f=u^{\mathrm{\prime }}+u\in L^2(0,1)$. But $u^{\mathrm{\prime }\mathrm{\prime }}(t)=-\frac{3}{16}{t}^{-5\mathrm{/}4}$ has nonintegrable square near zero, so $u\mathrm{\notin }{H}^2(0,1)$.

The exact equation integrated over $(0,\tau )$ is $\frac{u(\tau )-u(0)}{\tau }+\frac{1}{\tau }{\int }_0^{\tau }u(t)dt=\frac{1}{\tau }{\int }_0^{\tau }f(t)dt\mathrm{.}$ The first CN step from [F1] is $\frac{U^1-U^0}{\tau }+\frac{U^1+U^0}{2}=\frac{1}{\tau }{\int }_0^{\tau }f(t)dt$. Subtract, use $U^0=u(0)=0$, and set $e^1=U^1-u(\tau )$. This gives $e^1(\frac{1}{\tau }+\frac{1}{2})=\frac{1}{\tau }{\int }_0^{\tau }{t}^{3\mathrm{/}4}dt-\frac{{\tau }^{3\mathrm{/}4}}{2}=(\frac{4}{7}-\frac{1}{2}){\tau }^{3\mathrm{/}4}=\frac{1}{14}{\tau }^{3\mathrm{/}4}\mathrm{.}$

Solving the scalar equation in step 1.2 yields $e^1={\tau }^{7\mathrm{/}4}\mathrm{/}(14(1+\tau \mathrm{/}2))$. Along $\tau =1\mathrm{/}N\to 0$, $\mathrm{\mid }{e}^1\mathrm{\mid }\mathrm{/}{\tau }^2={\tau }^{-1\mathrm{/}4}\mathrm{/}(14(1+\tau \mathrm{/}2))\to \mathrm{\infty }$; the maximum nodal error is at least $\mathrm{\mid }{e}^1\mathrm{\mid }$, so it cannot be $O({\tau }^2)$.

The four triangles each have area $1\mathrm{/}4$, have disjoint interiors, and meet in whole edges or vertices; their union is the closed unit square. The piecewise barycentric coordinate $\varphi$ is zero on the boundary and agrees on shared edges. An affine polynomial on each triangle is determined by its three vertex values; matching traces make shared vertex values agree. Thus a global continuous P1 function is determined by the five mesh vertex values. The four boundary values vanish, leaving only the center value, so $V_h=\mathrm{span}\{\varphi \}$. By [F3], $\varphi \in H_0^1(\mathrm{\Omega })$; $u_0={\varphi }^2$ is also a continuous zero-boundary piecewise polynomial and lies in $H_0^1(\mathrm{\Omega })$. The norm and form hypotheses follow from [F4].

On a triangle $K$ of area $A=1\mathrm{/}4$, the affine map from the reference triangle has determinant magnitude $2A$. If $\lambda$ is the barycentric coordinate of the center vertex, direct integration gives ${\int }_K{\lambda }^{m}dx=2A{\int }_0^1{\int }_0^{1-x}(1-x-y{)}^{m}dydx=2A\mathrm{/}((m+1)(m+2))$ for $m=1,2,3$. Summing over the four triangles yields ${\int }_{\mathrm{\Omega }}\varphi =1\mathrm{/}3$, ${\int }_{\mathrm{\Omega }}{\varphi }^2=1\mathrm{/}6$, and ${\int }_{\mathrm{\Omega }}{\varphi }^3=1\mathrm{/}10$. The altitude from the center to each opposite boundary edge is $1\mathrm{/}2$, so the barycentric coordinate has gradient magnitude $2$ and $\mathrm{\mid }\mathrm{\nabla }\varphi {\mathrm{\mid }}^2=4$ on each triangle. Hence $a(\varphi ,\varphi )=4$ and $a({\varphi }^2,\varphi )=2{\int }_{\mathrm{\Omega }}\varphi \mathrm{\mid }\mathrm{\nabla }\varphi {\mathrm{\mid }}^2=8\mathrm{/}3$.

Write the orthogonal projection as $P_h{u}_0=p\varphi$. Its defining equation $(u_0-p\varphi ,\varphi {)}_{L^2}=0$ gives $p=(1\mathrm{/}10)\mathrm{/}(1\mathrm{/}6)=3\mathrm{/}5$. By [F5], the Ritz projection is $R_h{u}_0=r\varphi$ with $a(u_0-r\varphi ,\varphi )=0$, hence $r=(8\mathrm{/}3)\mathrm{/}4=2\mathrm{/}3$.

Therefore $P_h{u}_0-R_h{u}_0=(3\mathrm{/}5-2\mathrm{/}3)\varphi =-\varphi \mathrm{/}15\mathrm{\ne }0$, since $\varphi (c)=1$. If $u_h(0)=P_h{u}_0$, the error-splitting variable $\theta (0)=u_h(0)-R_h{u}_0$ is exactly this nonzero function.

The four triangles meet along whole edges and vertices and cover the closed square. On each triangle the center barycentric coordinate defines $\varphi$; these traces match, and its boundary trace is zero. An affine P1 function is determined by its three vertex values on each triangle, so the global zero-boundary P1 space has only the center value free and equals $\mathrm{span}\{\varphi \}$. Thus [F2] puts this trial space in $H_0^1(\mathrm{\Omega })$; [F3] supplies the Hilbert and form hypotheses. The function $u=q{\varphi }^2$ is also a zero-boundary conforming piecewise polynomial and is in the energy space. Since $q$ is smooth, $u^{\mathrm{\prime }}=q^{\mathrm{\prime }}{\varphi }^2\in V$; therefore $v\mapsto (u^{\mathrm{\prime }},v{)}_{L^2}$ and $v\mapsto a(u,v)$ are continuous functionals on $V$, and the defined load is square-integrable. By [F5] the initial Ritz value is defined.

On each triangle $K$, $\mathrm{\mid }K\mathrm{\mid }=1\mathrm{/}4$, the center barycentric coordinate $\lambda$ satisfies ${\int }_K{\lambda }^m=2\mathrm{\mid }K\mathrm{\mid }{\int }_0^1{\int }_0^{1-x}(1-x-y{)}^{m}dydx=2\mathrm{\mid }K\mathrm{\mid }\mathrm{/}((m+1)(m+2))$ for $m=1,2,3$ by affine change of variables and iterated integration. Hence ${\int }_{\mathrm{\Omega }}\varphi =1\mathrm{/}3$, ${\int }_{\mathrm{\Omega }}{\varphi }^2=1\mathrm{/}6$, and ${\int }_{\mathrm{\Omega }}{\varphi }^3=1\mathrm{/}10$. The altitude from the center to every boundary edge is $1\mathrm{/}2$, so $\mathrm{\mid }\mathrm{\nabla }\varphi {\mathrm{\mid }}^2=4$ on each triangle. Thus $a(\varphi ,\varphi )=4$ and $a({\varphi }^2,\varphi )=2{\sum }_K{\int }_K\varphi \mathrm{\mid }\mathrm{\nabla }\varphi {\mathrm{\mid }}^2=8\mathrm{/}3$. Also $({\varphi }^2,\varphi {)}_{L^2}=1\mathrm{/}10$.

By [F5], $R_{h}u(t)=(2\mathrm{/}3)q(t)\varphi$, since $a((2\mathrm{/}3)q\varphi ,\varphi )=q(8\mathrm{/}3)=a(q{\varphi }^2,\varphi )$. Therefore $\rho =q((2\mathrm{/}3)\varphi -{\varphi }^2)$ and $({\rho }^{\mathrm{\prime }},\varphi {)}_{L^2}=q^{\mathrm{\prime }}((2\mathrm{/}3)(1\mathrm{/}6)-1\mathrm{/}10)=q^{\mathrm{\prime }}\mathrm{/}90$.

Testing the semidiscrete equation with $\varphi$ and using step 1.2 gives $c^{\mathrm{\prime }}\mathrm{/}6+4c=q^{\mathrm{\prime }}\mathrm{/}10+8q\mathrm{/}3.$ For $\theta =(c-2q\mathrm{/}3)\varphi$, direct substitution yields $({\theta }^{\mathrm{\prime }},\varphi {)}_{L^2}+a(\theta ,\varphi )=c^{\mathrm{\prime }}\mathrm{/}6-q^{\mathrm{\prime }}\mathrm{/}9+4c-8q\mathrm{/}3=-q^{\mathrm{\prime }}\mathrm{/}90=-({\rho }^{\mathrm{\prime }},\varphi {)}_{L^2}\mathrm{.}$ This is [F6] for the sole basis vector, hence for every $v_h\in V_h$ by linearity. The calculation uses the initial Ritz projection, so $\theta (0)=0$.