Orthogonal cubic spline collocation method for the Cahn–Hilliard equation

Abstract

The Cahn–Hilliard equation plays an important role in the phase separation in a binary mixture. This is a fourth order nonlinear partial differential equation. In this paper, we study the behaviour of the solution by using orthogonal cubic spline collocation method and derive optimal order error estimates. We discuss some computational experiments by using monomial basis functions in the spatial direction and RADAU 5 time integrator. The method we present here is better in terms of stability, efficiency and conditioning of the resulting matrix. Since no integrals to be evaluated or approximated, it behaves better than finite element method.

Introduction

We consider the one spatial dimensional Cahn–Hilliard equation:

$$\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\gamma \frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^4}=\frac{{\mathrm{\partial }}^2\varphi (u)}{\mathrm{\partial }{x}^2}\text{,}(x\text{,}t)\in I\times (0\text{,}T]$$

with the initial condition

$$u(x\text{,}0)=u_0(x)\text{,}x\in I\text{,}$$

and boundary conditions

$$\begin{array}{cc}\hfill & u(0\text{,}t)=u(1\text{,}t)=0\text{,}t\in (0\text{,}T]\text{,}\hfill \\ \hfill & \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}(0\text{,}t)=\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}(1\text{,}t)=0\text{,}t\in (0\text{,}T]\text{,}\hfill \end{array}$$

where γ > 0, ϕ(u) = γ₂u³ + γ₁u² + γ₀u, γ₂ > 0, I = (0, 1) and 0 < T < ∞.

Eq. (1.1) arises in a variety of applications such as phase transition in material science. We refer the reader to [5] and the references there in. In this paper, we use second–order splitting procedure combined with orthogonal spline collocation method for Eq. (1.1) and derive optimal error estimates. Since this method is much superior to B-splines in terms of stability, efficiency and conditioning of the resulting matrix. Compared to finite element method (FEM) the calculation of the coefficients of the mass and stiffness matrices determining the approximate solution is very fast since no integrals need to be evaluated or approximated. We discuss numerical experiments using monomial basis functions and RADAU 5 time integrator.

Earlier, mixed methods in combination with orthogonal spline collocation methods used to fourth order evolution equations by Li et al. [9], Manickam et al. [12], [13], Danumjaya and Pani [3]. In the context of Cahn–Hilliard equation, Elliott et al. [6] discussed a second order splitting combined with lumped mass finite element method and derived optimal error estimates.

We split Eq. (1.1) by setting v = γu_xx − ϕ(u) then we obtain the following system:

$$u_t+v_{\mathit{xx}}=0\text{,}(x\text{,}t)\in I\times (0\text{,}T]\text{,}$$

$$\gamma u_{\mathit{xx}}-v-\varphi (u)=0\text{,}(x\text{,}t)\in I\times (0\text{,}T]$$

with initial condition

$$u(x\text{,}0)=u_0(x)\text{,}x\in I\text{,}$$

and the boundary conditions

$$\begin{array}{cc}\hfill & u(0\text{,}t)=u(1\text{,}t)=0\text{,}t\in (0\text{,}T]\text{,}\hfill \\ \hfill & v(0\text{,}t)=v(1\text{,}t)=0\text{,}t\in (0\text{,}T]\text{.}\hfill \end{array}$$

We use orthogonal cubic spline collocation method for the system (1.4), (1.5), (1.6), (1.7) in the spatial direction to compute the approximate solutions using monomial basis functions. Let $\{x_i{\}}_{i=1}^{N+1}$ denote a partition of $\overline{I}=[0\text{,}1]$ with

$$\begin{array}{cc}\hfill & 0=x_1<x_2<\cdots <x_{N+1}=1\text{,}\hfill \\ \hfill & I_j=(x_j\text{,}{x}_{j+1})\text{,}{h}_j=x_{j+1}-x_j\text{,}j=1\text{,}2\text{,}3\text{,}\dots \text{,}N\hfill \end{array}$$

and

$$h=\underset{1⩽j⩽N}{\max }{h}_j\text{.}$$

Assume that the partition is quasi-uniform, i.e., there exists a finite positive constant σ such that

$$\underset{1⩽j⩽N}{\max }\left(\frac{h}{h_j}\right)⩽\sigma \text{.}$$

We define a finite dimensional subspace ${\mathcal{H}}_3$ as

$${\mathcal{H}}_3=\{\chi \in C^1(\overline{I}):\chi {|}_{{\overline{I}}_j}\in P_3\text{,}j=1\text{,}2\text{,}\dots \text{,}N\mathrm{and}\chi (0)=\chi (1)=0\}\text{,}$$

where P₃ denotes the set of all cubic polynomials. Let $\{{\lambda }_k{\}}_{k=1}^2$ denote the roots of the Legendre polynomial of degree 2 i.e., $\left({\lambda }_1=\frac{1}{2}\left(1-\frac{1}{\sqrt{3}}\right)\text{,}{\lambda }_2=\frac{1}{2}\left(1+\frac{1}{\sqrt{3}}\right)\right)$. These are the nodes of the 2-point Gaussian quadrature rule on the interval I with corresponding weights w_k = 1/2, k = 1, 2. Now, we define the collocation point

$${\lambda }_{\mathit{jk}}=x_j+h_j{\lambda }_k\text{,}j=1\text{,}2\text{,}\dots \text{,}N\text{,}k=1\text{,}2\text{.}$$

Below, we define discrete innerproduct and its induced norm. For any $\varphi \text{,}\psi \in C^0(\overline{I})$, the discrete innerproduct is defined as

$$〈\varphi \text{,}\psi 〉=\sum _{j=1}^N〈\varphi \text{,}\psi {〉}_j\text{,}$$

where

$$〈\varphi \text{,}\psi {〉}_j=\frac{h_j}{2}\sum _{k=1}^2\varphi ({\lambda }_{\mathit{jk}})\psi ({\lambda }_{\mathit{jk}})\text{,}$$

and its induced discrete norm by

$$|\varphi {|}_D=〈\varphi \text{,}\varphi {〉}^{1/2}\text{.}$$

Lemma 1.1

For $w\text{,}z\in {\mathcal{H}}_3$,

$$-〈w^{″}\text{,}z〉=(w^{\prime }\text{,}{z}^{\prime })+\frac{1}{1080}\sum _{j=1}^N{w}_j^{(3)}{z}_j^{(3)}{h}_j^5=-〈z^{″}\text{,}w〉\text{,}$$

where $w_j^{(3)}$ (respectively, $z_j^{(3)}$), is the third derivative of w_j (respectively, z_j) which is constant on each subinterval I_j.

Note that when z = w with $w\in {\mathcal{H}}_3$, we have

$$\Vert w^{\prime }{\Vert }_{L^2}^2⩽-〈w^{″}\text{,}w〉\text{.}$$

For a proof of Lemma 1.1, we refer to Douglas and Dupont [4].

The outline of the paper is as follows. In Section 1, we introduce some notations and preliminaries. Section 2 deals with continuous-time orthogonal cubic spline collocation method for the solution of (1.4), (1.5), (1.6), (1.7). We establish optimal error estimates. Finally, Section 3 is devoted to numerical experiments. Here, we show that both theoretical order of convergence and numerically computed order of convergence are same.

Throughout this paper, C denotes a generic positive constant which is independent of the discretization parameter h which may have different values at different places.