Galerkin-Chebyshev Pseudo Spectral Method and a Split Step New Approach for a Class of Two dimensional Semi-linear Parabolic Equations of Second Order

Let us denote T_n(y) as the n^th degree Chebyshev polynomial and

PN=Span{T0(y),T1(y),…,TN(y)},VN={v∈PN:v(∓)=0}.$$\begin{array}{} \displaystyle P_{N}=Span\lbrace T_{0}(y),T_{1}(y),\ldots ,T_{N}(y)\rbrace , V_{N}=\lbrace v\in P_{N}:v(\mp )=0\rbrace . \end{array}$$(11)

Now, the standard Chebyshev-Galerkin approximation to equations (6) and (7) is nothing but finding X_N: [t_n,t_n+1] → VN2$\begin{array}{} V_{N}^{2} \end{array}$ and ϕ_N: [t_n,t_n+1] → VN2$\begin{array}{} V_{N}^{2} \end{array}$ such that

(X(tn+1)−X(tn),v)w−(G(tn+1)−G(tn),v)w=0,∀v∈VN2,tn≤t≤tn+1,$$\begin{array}{} \displaystyle (\mathrm{X}(t_{n+1})-\mathrm{X}(t_{n}),v)_{w}-(G(t_{n+1})-G(t_{n}),v)_{w}=0, \forall v\in V_{N}^{2}, t_{n}\le t\le t_{n+1}, \end{array}$$(12) and

(∂tϕN,v)w+k(t)(ΔϕN,v)w+(h(x,y,t,v)w=0,∀v∈VN2,tn≤t≤tn+1,(ϕN(x,y,t0)−ϕ0,v)w=0,∀v∈VN2,$$\begin{array}{} \displaystyle \begin{split} &{}(\partial _{t}\phi _{N},v)_{w}+k(t)(\Delta \phi_{N},v)_{w}+(h(x,y,t,v)_{w}=0, \forall v\in V_{N}^{2}, t_{n}\le t\le t_{n+1},\\ \displaystyle &{}\qquad \qquad \qquad (\phi _{N}(x,y,t_{0})-\phi _{0},v)_{w}=0, \forall v\in V_{N}^{2}, \end{split} \end{array}$$(13)

where inner product is defined in Lw2$\begin{array}{} \displaystyle L_{w}^{2} \end{array}$ (Ω) as (u,v)_w = ∫_Ωuvwdxdy, w = (1–x²)^–1/2(1– y²)^–1/2, ∥u∥w=∥u∥Lw2=(u,u)w1/2.$\begin{array}{} \Vert u\Vert _{w}=\Vert u\Vert _{L_{w}^{2}}=(u,u)_{w}^{1/2}. \end{array}$

We need the following lemma for our algorithm and proof of Lemma 1 can be found in Shen et al.[9]

Therefore, VN2$\begin{array}{} V_{N}^{2} \end{array}$ = {θ_k(x)θ_j(y), j,k = 0,1,…, N–2}. Let us assume that

uN(x,y,t)=∑k,j=0N−2Ckj(t)θk(x)θj(y).$$\begin{array}{} \displaystyle u^{N}(x,y,t)=\sum_{k,j=0}^{N-2}C_{kj}(t)\theta_{k}(x)\theta_{j}(y). \end{array}$$(17)

Then, for t∈ [t_n,t_n+1]

XN(x,y,t)=∑k,j=0N−2Dkj(t)θk(x)θj(y)withXN(x,y,tn)=u(x,y,tn),$$\begin{array}{} \displaystyle \mathrm{X}^{N}(x,y,t)=\sum_{k,j=0}^{N-2}{D_{kj}(t)\theta _{k}(x)\theta _{j}(y)} \, \mathrm{with}\, \mathrm{X}^{N}(x,y,t_{n})=u(x,y,t_{n}), \end{array} $$(18)

and

ϕN(x,y,t)=∑k,j=0N−2Ekj(t)θk(x)θj(y)withϕN(x,y,tn)=XN(x,y,tn+1).$$\begin{array}{} \displaystyle \phi^{N}(x,y,t)=\sum_{k,j=0}^{N-2}{E_{kj}(t)\theta_{k}(x)\theta_{j}(y)} \, \mathrm{with}\, \phi^{N}(x,y,t_{n})=\mathrm{X}^{N}(x,y,t_{n+1}). \end{array}$$(19)

Inserting (18), (19) into (12), (13) respectively, and using

v=θl(x)θs(y)l,s=0,1,2,…N−2.$$\begin{array}{} \displaystyle v=\theta _{l}(x)\theta _{s}(y) l,s=0,1,2,\ldots N-2. \end{array}$$

We get

(X(tn+1)−X(tn),θl(x)θs(y))w−(G(tn+1)−G(tn),θl(x)θs(y))w=0,tn≤t≤tn+1$$\begin{array}{} \displaystyle (\mathrm{X}(t_{n+1})-\mathrm{X}(t_{n}),\theta _{l}(x)\theta _{s}(y))_{w}-(G(t_{n+1})-G(t_{n}),\theta _{l}(x)\theta _{s}(y))_{w}=0, t_{n}\le t\le t_{n+1} \end{array}$$(20)

and

(∂tϕ^N,θl(x)θs(y))w+k(tn)(Δϕ^N,θl(x)θs(y))w=0,tn≤t≤tn+1(ϕ^N(x,y,t0)−ϕ^0,θl(x)θs(y))w=0.$$\begin{array}{} \displaystyle {}(\partial _{t}{\hat{\phi} }_{N},\theta _{l}(x)\theta _{s}(y))_{w}+k(t_{n})(\Delta {\hat{\phi} }_{N},\theta_{l}(x)\theta_{s}(y))_{w}=0, t_{n}\le t\le t_{n+1}\\ \displaystyle {}\qquad\qquad ({\hat{\phi} }_{N}(x,y,t_{0})-{\hat{\phi} }_{0},\theta_{l}(x)\theta_{s}(y))_{w}=0. \end{array}$$(21)

Then in the light of Lemma 1, we have

BDkj(tn+1)B−BDkj(tn)B−(G(tn+1)−G(tn))π=0,$$\begin{array}{} \displaystyle BD_{kj}(t_{n+1})B-BD_{kj}(t_{n})B-(G(t_{n+1})-G(t_{n}))\pi =0, \end{array}$$(22)

BddtEkj(tn+1)B−k(tn)(BEkj(tn+1)AT+AEkj(tn+1)B)=0,$$\begin{array}{} \displaystyle B\frac{d}{dt}E_{kj}(t_{n+1})B-k(t_{n})(BE_{kj}(t_{n+1})A^{T}+AE_{kj}(t_{n+1})B)=0, \end{array}$$(23)

BCkj(t0)B=U0,(U0)kj=(u0,θk(x)θj(y))w,$$\begin{array}{} \displaystyle BC_{kj}(t_{0})B=U_{0}, (U_{0})_{kj}=(u_{0},\theta _{k}(x)\theta _{j}(y))_{w}, \end{array}$$(24)

where A = (a_kj), B = (b_kj) k,j = 0,1,2,…, N–2.

We now present the following theorem for weighted Lw2$\begin{array}{} L_{w}^{2} \end{array}$ estimate for the discrete problem without proof. We refer to Liu et al. [5] for details.

where

Hws(Ω)={u:∥u∥s,w<∞},∥u∥s,w=∑|r|≤s∥Dru∥w212,|u|s,w=∑|r|=s∥Dru∥w212,$$\begin{array}{} \displaystyle H_{w}^{s}(\Omega )=\lbrace u:\Vert u\Vert _{s,w}<\infty \rbrace , \Vert u\Vert _{s,w}=\left(\sum_{\vert r\vert \le s}{\Vert D^{r}u\Vert _{w}^{2}}\right)^{\frac{1}{2}}, \vert u\vert _{s,w}=\left(\sum_{\vert r\vert =s}{\Vert D^{r}u\Vert _{w}^{2}}\right)^{\frac{1}{2}}, \end{array}$$(26)

here, D^ru = ∂|r|∂xr1∂yr2u,$\begin{array}{} \frac{\partial ^{\vert r\vert }}{\partial x^{r_{1}}\partial y^{r_{2}}}u, \end{array}$r = (r₁,r₂) and |r| = r₁ + r₂. Now, we define the Bochner space L^p((t₀,T); Hsw$\begin{array}{} H_{s}^{w} \end{array} $ (Ω )) endowed with the following norm

∥u∥Lp((t0,T);Hsw(Ω))=∫t0T∥u∥s,wpdt1/p,1≤p<∞esssupt∈(t0,T)∥u∥s,w,p=∞.$$\begin{array}{} \displaystyle \Vert u\Vert _{L^{p}((t_{0},T);H_{s}^{w}(\Omega ))}= \begin{cases} \left(\int_{t_{0}}^{T}{\Vert u\Vert _{s,w}^{p}dt}\right)^{1/p},& 1\le p<\infty \\ \displaystyle \mathrm{ess\, sup}_{t\in (t_{0},T)}\Vert u\Vert _{s,w},& p=\infty . \end{cases} \end{array} $$(27)

Let H0,w1={u∈Hw1:u(∂Ω)=0}.TheLw2(Ω)$\begin{array}{} H_{0,w}^{1}=\lbrace u\in H_{w}^{1}:u(\partial \Omega )=0\rbrace . {\rm The}\, L_{w}^{2}(\Omega ) \end{array}$ orthogonal projection PN0:H0,w1(Ω)→VN2$\begin{array}{} P_{N}^{0}:H_{0,w}^{1}(\Omega )\to V_{N}^{2} \end{array}$ is defined by

(PN0u−u,θ)w=0,∀θ∈VN2.$$\begin{array}{} \displaystyle (P_{N}^{0}u-u,\theta )_{w}=0, \forall \theta \in V_{N}^{2}. \end{array}$$(28)

We first briefly drive the model, the net electric charge density,ρ_e, in the EDL (electric double layer) can be related to the total potential field(Φ) by the Poisson equation

▽2Φ=−ρeϵ,$$\begin{array}{} \displaystyle \triangledown ^{2}\Phi =-\frac{\rho _{e}}{\epsilon }, \end{array} $$(29)

where ϵ is the permittivity of the solution. The total electric field Φ = Θ +Ψ is a linear superposition of the applied electric field generated by the electrodes, Θ (z), and the electric field due to the net charge distribution in the EDL, ψ(x,y), for details Movahed [10]. Consider the above equation with the standard electrokinetic model in Refs. [10,11]: Assume two dimensional, then RHS becomes a function of two variables,x and y. Substituting Φ = Θ + Ψ, we see that Eq.(18) can be separated very easily into the following two equations

▽2ψ=−ρeϵ,▽2Θ=0.$$\begin{array}{} \displaystyle \triangledown ^{2}\psi =-\frac{\rho _{e}}{\epsilon }, \triangledown ^{2}\Theta =0. \end{array} $$(30)

Therefore, the net electric charge density is independent of the external electric field and is determined by ψ. According to the quasi-equilibrium conditions, the ion density variation in the EDL obeys the Boltzmann distribution. So, for a symmetric electrolyte of equal valence (ℤ = ℤ⁺ = –ℤ^–), then the first part in Eq.(30) can be expressed explicitly by the PB equation as follows

▽2ψ=1ϵ2en0sinheZψkbTm,$$\begin{array}{} \displaystyle \triangledown ^{2}\psi =\frac{1}{\epsilon }\left[ 2en_{0}\sinh \left(\frac{e\mathbb{Z}\psi }{k_{b}T_{m}}\right)\right] , \end{array}$$(31)

where n₀ is the bulk number concentration of ions in the electrolyte solution, T_m is the absolute average temperature, k_b is the Boltzmann constant, and e is the elementary electronic charge. If use the dimensionless variables and parameters as in Refs. [10], [11], we obtain the dimensionless nonlinear Poisson-Boltzman. We further assume electrical field two-dimensional, in a two dimensional Cartesian coordinates system equation becomes

∂2ψ¯∂x2+∂2ψ¯∂y2=K2sinh⁡(ψ¯).$$\begin{array}{} \displaystyle \frac{\partial ^{2}\bar{\psi }}{\partial x^{2}}+\frac{\partial ^{2}\bar{\psi }}{\partial y^{2}}=K^{2}\sinh (\bar{\psi }). \end{array}$$(32)

This equation subject to the following boundary conditions (These conditions are widely suggested)

i)ψ¯(x,0)=ψ¯(y,0)=ψ¯(x,1)=ψ¯(y,1)=ψ0,ii)∂ψ¯∂x|x=0=0,∂ψ¯∂y|y=0=0,ψ¯(x,1)=ψ¯(y,1)=ψ0.$$\begin{array}{} \displaystyle {} i)\, \bar{ \psi }(x,0)=\bar{\psi }(y,0)=\bar{\psi }(x,1)=\bar{\psi}(y,1)=\psi _{0,}\\ \displaystyle {} ii)\, \frac{\partial \bar{\psi }}{\partial x}\vert _{x=0}=0,\frac{\partial \bar{\psi }}{\partial y}\vert _{y=0}=0,\bar{\psi}(x,1)=\bar{\psi }(y,1)=\psi _{0}. \end{array}$$(33)

Direct application of spectral method, requires solving highly nonlinear algebraic system which is not possible to solve numerically for large number of base elements. This motivates the development of pseudo-transient continuation approaches for solving the nonlinear PB equation (for details see Refs. [12, 13, 14, 15], the method authors used in this study are mainly finite difference and finite element methods. Basically, a pseudo-transient variation is introduced to convert Eq.(32) from the time-independent form to a time-dependent form

∂ψ¯∂t=∂2ψ¯∂x2+∂2ψ¯∂y2−K2sinh⁡(ψ¯).$$\begin{array}{} \displaystyle \frac{\partial \bar{\psi }}{\partial t}=\frac{\partial ^{2}\bar{\psi }}{\partial x^{2}}+\frac{\partial ^{2}\bar{\psi }}{\partial y^{2}}-K^{2}\sinh (\bar{\psi }). \end{array}$$(34)

We need now an initial solution, this could be the electrostatic potential obtained from a linearized PB equation or trivially ψ̄ = 0. We then numerically integrate Eq.(34) for a sufficiently large time period. The solution of the original nonlinear PB Eq.(32) can be recovered eventually by the steady state solution of the pseudo-time dependent process from Eq.(34).

However, there exists a real difficulty involved in the numerical integration of the time dependent NPB equation (34). Generally speaking, since a long time integration is required to reach the steady state for Eq.(34), explicit time stepping methods are usually not efficient for pseudo-transient continuation approaches, see for example [12, 13, 14, 15]. In the literature, semi-implicit time stepping methods [12, 13] are commonly used to solve time dependent NPB equation (34) so that a large time step could be used for a stable simulation. Nevertheless, a fully implicit time integration method has never been constructed for solving the classical nonlinear PB equation (34) by the method of Galerkin spectral method. Here we show that in Eqs.(12)-(13) and (10) that time splitting method enable us to use the implicit time integration. We first tested our algorithm for sinh (ψ̄) ≈ ψ̄ in which case Eq.(24) becomes linear and exact analytical solution is possible [6]. We found that L_w error less than 10^–8 for N = 10 and 0 < t < 10. Therefore, we expect similar error for nonlinear case. In Fig.1, we see that steady state case has been easily obtained. Furthermore, we are only able to show on figure first three approximations, because the other approximations are very small to show on the figure, which shows the power of our approximations.

Fig. 1

In Fig.2, we show the centerline electrostatic potential which shows the development of the profiles. We also show the three dimensional profile of electrostatic potential for t = 2.5. Figure 3 indicates the steady state solution, because any increase in the value of t, does not cause any change in the electrostatic potential.

Fig. 2

Fig. 3

Since the available experimental results for the Poisson-Boltzman equation are on [0,a], we obtained our results on this region. However, we note here that, in order to apply Chebyshev-Galerkin spectral method, we first transformed the given region to [–1,1], then we obtain the solution. Finally we applied inverse transformation to the obtained solution in order to express the solution in the original region. We also note that integration of equation in Eq. (12), we use Gauss-Chebyshev quadrature method.

We consider in this case k(t) = 1,a = c = 0, b = d = 1,α² = 1, f(u) = 11+u2$\begin{array}{} \displaystyle \frac{1}{1+u^{2}} \end{array}$ +h(x,y,t) and initial and boundary conditions

u(x,y,0)=sin⁡(πx)sin⁡(πy),onΩ=(0,1)×(0,1),$$\begin{array}{} \displaystyle u(x,y,0)=\sin (\pi x)\sin (\pi y), \, \mathrm{on}\, \Omega =(0,1)\times(0,1) , \end{array}$$(35)

u(−1,y,t)=u(1,y,t)=u(x,−1,t)=u(x,1,t)=0,fort>0$$\begin{array}{} \displaystyle u(-1,y,t)=u(1,y,t)=u(x,-1,t)=u(x,1,t)=0, \, \mathrm{for}\, t\gt0 \end{array}$$(36)

where the source function is chosen the such that exact solution is

u(x,y,t)=e−tsin⁡(πx)sin⁡(πy).$$\begin{array}{} \displaystyle u(x,y,t)=e^{-t}\sin (\pi x)\sin (\pi y). \end{array}$$(37)

In this case, according to our time splitting scheme, on [t_n, t_n+1] we get

∂χ∂t=1(1+χ2),withχ(x,y,tn)=u(x,y,tn),t∈[tn,tn+1],$$\begin{array}{} \displaystyle \frac{\partial \chi }{\partial t}=\frac{1}{(1+\chi ^{2})},\, \mathrm{with}\, \chi (x,y,t_{n})=u(x,y,t_{n}), t\in \lbrack t_{n},t_{n+1}\rbrack , \end{array}$$(38)

∂ϕ∂t=(∂2ϕ∂x2+∂2ϕ∂y2)+h(x,y,t),withϕ(x,y,tn)=χ(x,y,tn+1),t∈[tn,tn+1],$$\begin{array}{} \displaystyle \frac{\partial \phi }{\partial t}=(\frac{\partial ^{2}\phi }{\partial x^{2}}+\frac{\partial ^{2}\phi }{\partial y^{2}})+h(x,y,t),\, \mathrm{with}\, \phi (x,y,t_{n})=\chi (x,y,t_{n+1}), t\in \lbrack t_{n},t_{n+1}\rbrack , \end{array}$$(39)

where,

h(x,y,t)=−e−3t(sin⁡(πx)sin⁡(πy))3(1−2π2)−e−tsin⁡(πx)sin⁡(πy)(1−2π2)−1e−2t(sin⁡(πx)sin⁡(πy))2+1.$$\begin{array}{} \displaystyle h(x,y,t)=\dfrac{-e^{-3t}(\sin (\pi x)\sin (\pi y))^{3}(1-2\pi ^{2})-e^{-t}\sin (\pi x)\sin (\pi y)(1-2\pi ^{2})-1}{e^{-2t}(\sin (\pi x)\sin (\pi y))^{2}+1}. \end{array}$$

Integration of equation (38) and using the Cardana formula, we get

χ(x,y,tn+1)−χ(x,y,tn)=12(12tn+1+44+9(tn+1)2)13−2(12tn+1+44+9(tn+1)2)13−12(12tn+44+9(tn)2)13−2(12tn+44+9(tn)2)13.$$\begin{array}{} \displaystyle &{}\chi (x,y,t_{n+1})-\chi (x,y,t_{n})\\ \displaystyle &{}\qquad=\frac{1}{2}(12t_{n+1}+4\sqrt{4+9(t_{n+1})^{2}})^{\frac{1}{3}}-\frac{2}{(12t_{n+1}+4\sqrt{4+9(t_{n+1})^{2}})^{\frac{1}{3}}}\\ \displaystyle &{}\qquad-\frac{1}{2}(12t_{n}+4\sqrt{4+9(t_{n})^{2}})^{\frac{1}{3}}-\frac{2}{(12t_{n}+4\sqrt{4+9(t_{n})^{2}})^{\frac{1}{3}}} . \end{array}$$(40)

Since we have solved nonlinear part analytically, we can apply implicit scheme for temporal discretization. Now, substituting (18) and (19) into (30) and (31), and applying the Chebyshev Galerkin method, we obtain simple expression temporal discretization which can be solved by implicit method as described earlier. Hence, we see that our method much more efficient than the usual spectral Galerkin method. Fig. 4 shows the approximation for the h(x,y,t) in terms of our base function in Eq.(14) for given fixed time. We see that the error is less than 10^–4 except for the initial and end regions where h(x,y,t) is not zero; whereas our base function satisfies the boundary condition zero at the end-point. Fig. 5 shows the L_w error for different N.

Fig. 4

Fig. 5