# How to verify correctness of my numerical method for Allen-Cahn equation?

$$\label{Parabolic} \frac{\partial \phi(\mathbf{x},t)}{\partial t} - \Delta\phi(\mathbf{x},t)+\frac{f(\phi(\mathbf{x},t))}{\epsilon^2}=0$$

$$\label{boundary} \frac{\partial \phi(\mathbf{x},t)}{\partial \mathbf{n}}=0, \ \text{on} \ \partial \Omega \times [0,T]$$

$$\label{initail} \phi(\mathbf{x},0) = \phi_{0}(\mathbf{x}) \ \text{in} \ \ \Omega$$ Where $f(\phi)=\phi-\phi^{3}$, $\epsilon$ is a positive constant, and $\mathbf{n}$ is the outward normal vector at the domain boundary.

I use finite element method to solve this equation:

Firstly, I discrete the equation in time space $$\label{timediscrete} \frac{\phi^{n+1}-\phi^{n}}{\tau}-\Delta\phi^{n+1}=-\frac{f(\phi^{n})}{\epsilon^2}$$ This is not a strict back Euler method, we get a elliptic equation with homogeneous neumann boundary condition： $$\label{ellipse} \frac{\phi^{n+1}}{\tau}-\Delta\phi^{n+1}=-\frac{f(\phi^{n})}{\epsilon^2}+\frac{\phi^{n}}{\tau}$$

Suppose $V_{h}\in H^{1}(\Omega)$ is the finite element space, the problem can be stated as follows: find $\phi_{h}(\cdot,t)\in V_{h}$ s.t. $$\label{variation} \frac{1}{\tau}(\phi_h^{n+1},v_h)-a(\phi_h^{n+1},v_h)=(g^n,v_h)$$ where $v_h \in V_h, a(u,v)=\int_{\Omega}\nabla u\cdot \nabla v dxdy,\ g^n=-\frac{f(\phi^{n})}{\epsilon^2}+\frac{\phi^{n}}{\tau}$.

Assume $$\phi_{h}(\mathbf{x},t)=\sum_{i=1}^{N}\mu_{i}(t)\varphi_{i}(\mathbf{x})$$, put this in the above and set $v_h=\varphi_{i}$, we can get the following equations:

$$\label{equationt} \sum_{i=1}^{N}(\frac{1}{\tau}(\varphi_i,\varphi_j)+a(\varphi_i,\varphi_j))\mu_i(t_{n})=(g^n,\varphi_j)$$ with initial condition: $$\label{equation0} \sum_{i=1}^{N}(\varphi_i,\varphi_j)\mu_i(0)=(\varphi_0,\varphi_j)$$ $j=1,2,\cdots,N$.

As I don't know the exact solution of the problem, How can I verify my numerical method? I use $k$th degree polynomial as the basis function, is there any theoretical convergence order of my numerical method?

• I would advise against simultaneously linearizing and discretizing in time. If you first e.g. discretize in time and then linearize using Newton's method then the method can be studied step-by-step and it will be easier to find theoretical error bounds. Regarding your current method: If you cannot find any exact solution then I'm afraid you are out of luck. – knl Dec 13 '15 at 10:59
• I am also worried about my numerical method. But I don't understand what linearize using Newton's method mean here. Can you please tell me more detail? Thank You! – yang Dec 13 '15 at 12:35
• First applying implicit Euler in time gives $$\tau^{-1} \phi^{n+1} - \Delta \phi^{n+1}+\frac{f(\phi^{n+1})}{\epsilon^2}=\tau^{-1} \phi^n.$$ This is now a nonlinear equation in $\phi^{n+1}$ and can be linearized by writing $\phi^{n+1} = \phi^{n+1}_0 + \tilde{\phi}^{n+1}$ where $\phi^{n+1}_0$ is some initial guess and $\tilde{\phi}^{n+1}$ is the new unknown. This gives $$\tau^{-1}\tilde{\phi}^{n+1}-\Delta \tilde{\phi}^{n+1}+\epsilon^{-2}(\tilde{\phi}^{n+1}-3(\phi_0^{n+1})^2 \tilde{\phi}^{n+1})=\tau^{-1}(\phi^n-\phi^{n+1}_0)+\Delta\phi^{n+1}_0-\epsilon^{-2}(\phi_0^{n+1}+(\phi_0^{n+1})^3)$$ – knl Dec 13 '15 at 15:56
• Thank you for answering me. There may be a little mistake in the last equation. In the last term $-\epsilon^{-2}(\phi_0^{n+1}+(\phi_0^{n+1})^3)$ may be $-\epsilon^{-2}(\phi_0^{n+1}-(\phi_0^{n+1})^3)$. But I wonder why you write $\phi^{n+1}=\phi_0^{n+1}+\tilde{\phi}^{n+1}$. Can only expand $f(\phi^{n+1})$ at $\phi_0^{n+1}$? that's to say: $f(\phi^{n+1})=f(\phi_0^{n+1})+f'(\phi_0^{n+1})(\phi^{n+1}-\phi_0^{n+1})$, We get $\tau^{-1}\phi^{n+1}-\Delta \phi^{n+1}+\epsilon^{-2}(\phi^{n+1}-3(\phi_0^{n+1})^2\phi^{n+1})=\tau^{-1}\phi^{n}-2\epsilon^{-2}(\phi_0^{n+1})^3$ – yang Dec 14 '15 at 1:19
• You are correct about the mistake. The idea of Newton's method is to repeatedly compute the perturbation $\tilde{\phi}^{n+1}$ and sum it to the initial guess until you get convergence. Thus, here you would have two iterations, one for time and (at each time step) one for the nonlinearity. This way you can verify both steps individually and since both methods (implicit Euler and Newton's) are well studied in the literature you can find some convergence rates and related theory for verification. – knl Dec 14 '15 at 13:54

1. Instead solving $Lu=0$ solve $Lu=f$ subject to inhomogeneous boundary condition $Bu=g$. Next, choose an arbitrary function $u$ and calculate $f=Lu$ and $g=Bu$ analytically.
For example for exact solution $u=e^{at}e^{bx+dy}$ on rectangular domain you use $$f = \left(a - b^2 - d^2 +\frac{1-e^{2at}e^{2(bx+dy)}}{\epsilon^2}\right)e^{at}e^{bx+dy}$$ subject to $$\phi_x(x_0,y,t)=be^{at}e^{bx_0+dy},\\ \phi_x(x_n,y,t)=be^{at}e^{bx_n+dy},\\ \phi_y(x,y_0,t)=de^{at}e^{bx+dy_0},\\ \phi_y(x,y_n,t)=de^{at}e^{bx+dy_n},\\$$ and $$\phi(x,y,0) = e^{bx+dy} \ \text{in} \ \ \Omega$$
2. Instead of testing against an unknown exact solution you refine the grid and measure sort of Cauchy convergence, i.e. $e=\|u_h-u_{h/2}\|$ where $u_h$ is a solution obtained on grid of size\step $h$. Measure the rate of convergence of this error which will be the same as for $e=\|u-u_h\|$ where $u$ is exact solution if known. You will need odd number of points of the grid to be able to match between subsequent grids (the finer one has twice more points, but you test only for those that a common for both grids).