Analytic approximation of solutions of parabolic partial differential equations with variable coefficients

TL;DR

This paper develops an explicit solution method for one-dimensional reaction-diffusion equations with variable coefficients, enabling accurate numerical solutions and extending classical heat polynomial techniques through transmutation operators.

Contribution

It introduces a new family of solutions derived via transmutation operators, facilitating explicit solutions and numerical approximations for variable-coefficient parabolic PDEs.

Findings

Solutions are explicitly constructed using transmutation operators.

The method achieves high accuracy in numerical solutions.

Applicable to various boundary conditions beyond Dirichlet.

Abstract

A complete family of solutions for the one-dimensional reaction-diffusion equation \[ u_{xx}(x,t)-q(x)u(x,t) = u_t(x,t) \] with a coefficient $q$ depending on $x$ is constructed. The solutions represent the images of the heat polynomials under the action of a transmutation operator. Their use allows one to obtain an explicit solution of the noncharacteristic Cauchy problem for the considered equation with sufficiently regular Cauchy data as well as to solve numerically initial boundary value problems. In the paper the Dirichlet boundary conditions are considered however the proposed method can be easily extended onto other standard boundary conditions. The proposed numerical method is shown to reveal good accuracy.

Tables

Table 1. Table 1: Maximal absolute and relative errors of the approximate solution and condition number of the matrix in ( 5.5 ) for the problem ( 6.1 )–( 6.3 ) obtained for different values of N 𝑁 N in ( 5.4 ).

$N$	Max. absolute error	Max. relative error	Cond. number
5	$2.3\cdot {10}^{-2}$	$6.2\cdot {10}^{-2}$	$55.7$
10	$1.4\cdot {10}^{-4}$	$5.5\cdot {10}^{-4}$	$1.82\cdot {10}^5$
15	$9.6\cdot {10}^{-8}$	$4.0\cdot {10}^{-7}$	$3.65\cdot {10}^9$
20	$2.0\cdot {10}^{-10}$	$8.6\cdot {10}^{-10}$	$1.76\cdot {10}^{14}$
23	$7.6\cdot {10}^{-13}$	$3.2\cdot {10}^{-12}$	$1.67\cdot {10}^{17}$
26	$1.8\cdot {10}^{-13}$	$7.8\cdot {10}^{-13}$	$2.59\cdot {10}^{23}$
29	$2.5\cdot {10}^{-12}$	$1.1\cdot {10}^{-11}$	$1.23\cdot {10}^{23}$
34	$1.7\cdot {10}^{-10}$	$7.3\cdot {10}^{-10}$	$1.37\cdot {10}^{25}$
39	$2.3\cdot {10}^{-9}$	$9.8\cdot {10}^{-9}$	$2.86\cdot {10}^{29}$
50	$4.7\cdot {10}^{-10}$	$2.1\cdot {10}^{-9}$	$1.25\cdot {10}^{41}$
75	$6.1\cdot {10}^{-11}$	$2.7\cdot {10}^{-10}$	$2.89\cdot {10}^{73}$
100	$2.8\cdot {10}^{-10}$	$1.2\cdot {10}^{-9}$	$5.63\cdot {10}^{105}$

Equations

$$u_{xx}(x,t)-q(x)u(x,t)=u_t(x,t)$$

$$f^{\prime \prime }-q(x)f=0$$

$$X^{(0)}\equiv 1,X^{(n)}(x)$$

$${\tilde{X}}^{(0)}\equiv 1,{\tilde{X}}^{(n)}(x)$$

$${\phi }_k(x)=\{\begin{array}{ll}f(x)X^{(k)}(x),& k\text{ odd}\\ f(x){\tilde{X}}^{(k)}(x),& k\text{ even}\end{array}$$

$$ATv=TBv$$

$$Tv(x):=v(x)+{\int }_{-x}^{x}K(x,s)v(s)ds.$$

$$K_{xx}(x,s)-q(x)K(x,s)=K_{ss}(x,s),$$

$$K(x,x)=\frac{\alpha }{2}+\frac{1}{2}{\int }_0^{x}q(y)dy,x\in [-b,b],$$

$$K(x,-x)=\frac{\alpha }{2},x\in [-b,b].$$

$$H_{uv}(u,v)=q(u+v)H(u,v),$$

$$H(u,0)=\frac{\alpha }{2}+\frac{1}{2}{\int }_0^{u}q(s)ds,H(0,v)=\frac{\alpha }{2}.$$

$$T[x^k]={\phi }_k(x),\forall k\in {\mathbb{N}}_0.$$

$$v(0)=u(0),v^{\prime }(0)=u^{\prime }(0)-\alpha u(0),$$

$$\left|v^{(j)}(t)\right|\le C_1{\gamma }_1^{-2j}(2j)!,j=0,1,\dots$$

$$u_{xx}(x,t)-q(x)u(x,t)$$

$$u(0,t)$$

$$u_x(0,t)$$

$$\sum _{j=0}^{\infty }\left[\frac{F^{(j)}(t)}{(2j)!}\left({\phi }_{2j}(x)-\frac{\alpha }{2j+1}{\phi }_{2j+1}(x)\right)+\frac{G^{(j)}(t)}{(2j+1)!}{\phi }_{2j+1}(x)\right]$$

$$h_{xx}$$

$$h(0,t)$$

$$h_x(0,t)$$

$$h(x,t)=\sum _{k=0}^{\infty }\left(\frac{F^{(k)}(t)}{(2k)!}{x}^{2k}+\frac{G^{(k)}(t)-\alpha F^{(k)}(t)}{(2k+1)!}{x}^{2k+1}\right).$$

$$u(x,t)=Th(x,t)=\sum _{k=0}^{\infty }\left[\frac{F^{(k)}(t)}{(2k)!}\left({\phi }_{2k}(x)-\frac{\alpha }{2k+1}{\phi }_{2k+1}(x)\right)+\frac{G^{(k)}(t)}{(2k+1)!}{\phi }_{2k+1}(x)\right].$$

$$h_n(x,t)=n!\sum _{k=0}^{\left[\frac{n}{2}\right]}\frac{t^k{x}^{n-2k}}{k!(n-2k)!},n\in {\mathbb{N}}_0.$$

$$u_n(x,t)=n!\sum _{k=0}^{\left[\frac{n}{2}\right]}\frac{t^k{\phi }_{n-2k}(x)}{k!(n-2k)!},n\in {\mathbb{N}}_0$$

$$\frac{{\partial }^2{u}_n}{\partial x^2}-q(x)u_n(x,t)-\frac{\partial u_n}{\partial t}(x,t)=T\left(\frac{{\partial }^2{h}_n}{\partial x^2}-\frac{\partial h_n}{\partial t}(x,t)\right)=0.$$

$$\underset{\bar{\Omega }}{\max }\left|u(x,t)-\sum _{n=0}^N{a}_n{u}_n(x,t)\right|<\epsilon .$$

$$\underset{\bar{\Omega }}{\max }\left|h(x,t)-\sum _{n=0}^N{a}_n{h}_n(x,t)\right|<{\epsilon }_1.$$

$$\underset{\bar{\Omega }}{\max }\left|u(x,t)-\sum _{n=0}^N{a}_n{u}_n(x,t)\right|=\underset{\bar{\Omega }}{\max }\left|Th(x,t)-\sum _{n=0}^N{a}_{n}T{h}_n(x,t)\right|\le C{\epsilon }_1$$

$$u_{xx}(x,t)-q(x)u(x,t)=u_t(x,t),(x,t)\in \Omega$$

$$u(-b,t)={\psi }_1(t),u(b,t)={\psi }_2(t),t\in [0,\tau ]$$

$$u(x,0)=\phi (x),x\in [-b,b]$$

Full text

[TABLE]

with a coefficient $q$ depending on $x$ is constructed. The solutions represent the images of the heat polynomials under the action of a transmutation operator. Their use allows one to obtain an explicit solution of the noncharacteristic Cauchy problem for equation (0.1) with sufficiently regular Cauchy data as well as to solve numerically initial boundary value problems. In the paper the Dirichlet boundary conditions are considered however the proposed method can be easily extended onto other standard boundary conditions. The proposed numerical method is shown to reveal good accuracy.

1 Introduction

In the present work a complete system of solutions of a one-dimensional reaction-diffusion equation with a variable coefficient

[TABLE]

considered on $\Omega:=(-b,b)\times(0,\tau)$ is obtained. We assume that the potential $q\in C[-b,b]$ may be complex valued. The completeness of the system is with respect to the uniform norm in the closed rectangle $\bar{\Omega}$. The system of solutions is shown to be useful for uniform approximation of solutions of initial boundary value problems for (1.1)111In the paper the Dirichlet boundary conditions are considered however the proposed method can be easily extended onto other standard boundary conditions. as well as for explicit solution of the noncharacteristic Cauchy problem (see [2]) for (1.1) in terms of the formal powers arising in the spectral parameter power series (SPPS) method (see [8], [10], [13]).

The complete system of solutions is constructed with the aid of the transmutation operators relating (1.1) with the heat equation (see e.g., [4], [16], [5]). The possibility to construct complete systems of solutions by means of transmutation operators was proposed and explored in [4], and the approach developed in [4] requires the knowledge of the transmutation operators. In the present work using a mapping property of the transmutation operators discovered in [1] we show that the construction of the complete systems of solutions for equations of the form (1.1), representing transmuted heat polynomials, can be realized with no previous construction of the transmutation operator. Moreover, the use of the mapping property leads to an explicit solution of the noncharacteristic Cauchy problem for (1.1) with Cauchy data belonging to a Holmgren class [2].

We illustrate the implementation of the complete system of the transmuted heat polynomials by a numerical solution of an initial boundary value problem for (1.1). The approximate solution is sought in the form of a linear combination of the transmuted heat polynomials and the initial and boundary conditions are satisfied by a collocation method. A remarkable accuracy is achieved in few seconds using Matlab 2012 on a usual PC.

Besides this Introduction the paper contains five sections. In Section 2 we recall the transmutation operators and some of their properties. In Section 3 an explicit solution of the noncharacteristic Cauchy problem for (1.1) is obtained. In Section 4 we prove the completeness of the transmuted heat polynomials. In Section 5 the numerical method for solving initial boundary value problems for (1.1) implementing the transmuted heat polynomials is discussed. Section 6 presents a numerical illustration.

2 Transmutation operators and formal powers

2.1 System of recurrent integrals

Let $q$ be a continuous complex valued function defined on the segment $[-b,b]$. Throughout the paper we suppose that $f$ is a nonvanishing solution of the equation

[TABLE]

on $(-b,b)$ such that $f(0)=1$ and $f^{\prime}(0)=\alpha$ where $\alpha$ is a complex number. In [10] the existence of such solution was proved.

Consider two sequences of recurrent integrals (see [8], [10], [9])

[TABLE]

Definition 2.1**.**

The family of functions $\{\varphi_{k}\}_{k=0}^{\infty}$ constructed according to the rule

[TABLE]

is called the system of formal powers associated with $f$.

The formal powers arise in the spectral parameter power series (SPPS) representation for solutions of the one-dimensional Schrödinger equation (see [8], [10]).

2.2 The transmutation operator

For any $q\in C[-b,b]$ it is a well known result [16, Chapter 1] that there exists a function $(x,s)\mapsto K(x,s)$ defined on the domain $0\leq\left|s\right|\leq|x|\leq b$, continuously differentiable, such that the equality

[TABLE]

is valid for all $v\in C^{2}[-b,b]$, where $A:=\frac{\partial^{2}}{\partial x^{2}}-q$, $B:=\frac{\partial^{2}}{\partial x^{2}}$ and $T$ has the form of a second kind Volterra integral operator

[TABLE]

The operator $T$ is called transmutation operator. Moreover, the function $K$ is not unique and can be chosen so that $T[1]=f$ (see, e.g., [14]). When $q\in C^{1}[-b,b]$ such function $K$ is the unique solution of the Goursat problem

[TABLE]

For any $q\in C[-b,b]$ the kernel $K$ can be defined as K(x,s)=H\big{(}\frac{x+s}{2},\frac{x-s}{2}\big{)}, $|s|\leq|x|\leq b$, with $H$ being the unique solution of the Goursat problem

[TABLE]

If the potential $q$ is $n$ times continuously differentiable on $(-b,b)$, the kernel $K(x,s)$ is $n+1$ times continuously differentiable with respect to both independent variables.

The following mapping property of the operator $T$ is used throughout the paper.

Proposition 2.2** ([1]).**

[TABLE]

The inverse operator $T^{-1}$ also has the form of a second kind Volterra integral operator and satisfy the following correspondence of the initial values, see [11]

[TABLE]

where $v:=T^{-1}u$.

3 The noncharacteristic Cauchy problem for (1.1)

In this section an explicit solution of the noncharacteristic Cauchy problem for (1.1) in terms of the formal powers $\varphi_{k}$ is obtained. This fact is a direct consequence of the mapping property (2.3).

Definition 3.1** ([2]).**

For the positive constants $\gamma_{1}$, $\gamma_{2}$ and $C_{1}$, the Holmgren class $H(\gamma_{1},\gamma_{2},C_{1},t_{0})$ is the set of infinitely differentiable functions $v$ defined on $\left|t-t_{0}\right|<\gamma_{2}$ that satisfy

[TABLE]

for all $t\in\left|t-t_{0}\right|<\gamma_{2}$.

Proposition 3.2**.**

Let $q\in C[-b,b]$ and $u(x,t)$ be a solution of the noncharacteristic Cauchy problem

[TABLE]

where $F,G\in H(b,\tau,C,0)$, $C>0$. Then the series

[TABLE]

converges uniformly and absolutely for $\left|x\right|\leq r<b$ to the solution $u(x,t)$ where $\varphi_{k}$ are the formal powers (2.2).

Remark 3.3*.*

Note that the functions $\varphi_{2j}(x)-\frac{\alpha}{2j+1}\varphi_{2j+1}(x)$ coincide with the formal powers (2.2) constructed starting with the particular solution $g$ of (2.1) satisfying the initial conditions $g(0)=1$ and $g^{\prime}(0)=0$, see [12, Proposition 4.7].

Proof. Let $u(x,t)$ be a solution of ($\ref{eq cauchy})$. Consider the function $h:=T^{-1}u$, where the operator $T^{-1}$ is applied with respect to the variable $x$. The function $h$ is a solution of the heat equation, c.f., [4, Theorem 2.1.2] and due to (2.4) satisfies the following noncharacteristic Cauchy problem

[TABLE]

Since $G-\alpha F\in H(b,\tau,(1+\alpha)C,0)$, the solution of this problem is given by the absolutely and uniformly convergent series for $\left|x\right|\leq r<b$ (see, e.g., [2])

[TABLE]

Due to (2.3) we obtain

[TABLE]

This series converges uniformly and absolutely for $\left|x\right|\leq r<b$ due to the uniform boundedness of $T$ and of its inverse.

4 Transmuted heat polynomials

In this section a complete system of solutions of (1.1) is presented.

Consider the heat polynomials (see, e.g., [18]) defined by

[TABLE]

Due to (2.3) we obtain that the functions

[TABLE]

are solutions of (1.1) for all $-b<x<b$ and $t>0$. Indeed, we have that $u_{n}=Th_{n},$ $n\in\mathbb{N}_{0}$ and

[TABLE]

The completeness of the system of the heat polynomials with respect to the maximum norm proved in [4] and the uniform boundedness of $T$ and $T^{-1}$ imply the completeness of (4.2) in the space of classical solutions of (1.1). Thus, the following statement is true.

Theorem 4.1**.**

Let $u(x,t)$ be continuous in $\bar{\Omega}$ and satisfy (1.1) in $\Omega$. Then given $\varepsilon>0$ there exists $N\in\mathbb{N}$ and constants $a_{0},a_{1},\ldots,a_{N}$ such that

[TABLE]

Proof. Choose $\varepsilon>0$. Consider $h(x,t)=T^{-1}u(x,t)$. Due to the completeness of (4.1), for any $\varepsilon_{1}>0$ there exists $N\in\mathbb{N}$ and constants $a_{0},a_{1},\ldots,a_{N}$ such that

[TABLE]

Then

[TABLE]

where the constant $C$ is the uniform norm of $T$. The choice of $\varepsilon_{1}=\varepsilon/C$ finishes the proof.

5 Solution of initial boundary value problems for (1.1)

Consider the problem to find the solution of the equation

[TABLE]

subject to the Dirichlet boundary conditions

[TABLE]

and the initial condition

[TABLE]

where $\psi_{1}$, $\psi_{2}$ and $\varphi$ are continuously differentiable functions satisfying the compatibility conditions

[TABLE]

The problem (5.1)–(5.3) possesses a unique solution and depends continuously on the data (see, e.g., [17]).

The result of Theorem 4.1 suggests the following simple method to approximate the solution of problem (5.1)–(5.3). The approximate solution $\tilde{u}$ is sought in the form

[TABLE]

Since every $u_{n}$ is a solution of (5.1), their linear combination satisfies (5.1) as well. The coefficients $\left\{a_{n}\right\}_{n=0}^{N}$ are sought in such way that $\tilde{u}$ satisfy the initial and the boundary conditions approximately. For this we used the collocation method. $M$ points $\{(x_{i},t_{i})\}_{i=1}^{M}$ are chosen on the parabolic boundary $(\{-b\}\times[0,\tau])\cup([-b,b]\times\{0\})\cup(\{b\}\times[0,\tau])$ in order to construct a linear system of equations for the coefficients $\left\{a_{n}\right\}_{n=0}^{N}$

[TABLE]

The system (5.5) is the result of imposing the conditions (5.2), (5.3) onto the approximate solution (5.4). Using the pseudoinverse matrix the system (5.5) is solved, and the approximate solution (5.4) is computed on $\Omega$ using the obtained coefficients $\left\{a_{n}\right\}_{n=0}^{N}$ and the definition of the transmuted heat polynomials (4.2).

Needless to add that the same approach is applicable to other kinds of boundary conditions.

6 Numerical illustration

We present a numerical example of the application of the method described in the previous section. It reveals a remarkable accuracy with very little computational efforts. The implementation was realized in Matlab 2012.

On the first step a nonvanishing solution $f$ of (2.1) was computed using the SPPS method (see [10], [13]). The formal powers $\varphi_{k}$ were constructed like in [5] using the spapi and fnint Matlab routines from the spline Toolbox. Then the transmuted heat polynomials (4.2) were calculated. In order to obtain a unique solution of the system (5.5) $N+1$ equally spaced points on the parabolic boundary were chosen. Finally, the approximate solution (5.4) was computed on a mesh of $200\times 100$ points in the interior of the rectangle and compared with the corresponding exact solution.

Example 6.1*.*

Consider the initial Dirichlet problem

[TABLE]

The exact solution of this problem has the form

[TABLE]

The distribution of the absolute error of the approximate solution for $N=26$ is presented on Figure 1. The maximum absolute error of the approximate solution is of order $10^{-13}$.

It is often stated that boundary collocation methods (in particular, the heat polynomials method) lead to ill-conditioned systems of linear equations, see [3], [6], [15]. It is also the case for the proposed method. As is illustrated in Table 1, the condition number of the matrix in (5.5) grows rather fast. Nevertheless, the straightforward implementation of the proposed method presented no numerical difficulties. The convergence and the robustness of the method are illustrated in Table 1 where the maximum absolute and the maximum relative error of the approximate solution for different values of $N$ used for approximation (5.4) are presented. As one can appreciate, the convergence rate for small values of $N$ is exponential. And even taking values of $N$ much larger than the optimum one do not lead to any problem for collocation method nor to significant precision lost. Moreover, a simple test based on the accuracy of fulfilment of the initial and boundary conditions (6.2)–(6.3) can be utilized to estimate both the optimal $N$ and the accuracy of the obtained approximate solution.

Conclusions

A complete system of solutions of equation (1.1) is obtained. The solutions represent the images of the heat polynomials under the action of the transmutation operator. They are shown to be convenient for uniform approximation of solutions of initial boundary value problems for (1.1) as well as for explicit solution of the noncharacteristic Cauchy problem. Besides the Dirichlet boundary conditions considered in this paper the method is applicable to other standard boundary conditions. The complete system of solutions obtained can be used for solving moving and free boundary problems [7].

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

Research was supported by CONACYT, Mexico via the project 222478.