Extension of the operational Tau method for solving 1-D nonlinear transient heat conduction equations

Lemma 2.1 Ortiz and Samara, 1981

If $y_N(x)={\underline{a}}_N{X}_x$ , where a_N = (a₀,a₁, … ,a_N,0,…) and X_x = (1,x, … ,x^N,…)^T, then

1. xy_N(x) = a_NμX_x;

2. $\frac{d}{\mathit{dx}}{y}_N(x)={\underline{a}}_N\eta X_x$ .

Corollary 2.2 Ortiz and Samara, 1981

Generally, under assumptions of Lemma 2.1, we have

1. xⁱX_x = μⁱX_x;

2. $\frac{d^r}{{\mathit{dx}}^r}{X}_x={\eta }^r{X}_x$ .

If u_N(x,t) = u ϕ_x,t, then

1. ${\mathit{xu}}_N(x\text{,}t)=\underline{u}{\hat{\mu }}_x{\underline{\varphi }}_{x\text{,}t}$ , where ${\hat{\mu }}_x=\Phi {\tilde{\mu }}_x{\Phi }^{-1}\text{,}{\tilde{\mu }}_x=\mu \otimes I$ with the elements ${({\tilde{\mu }}_x)}_{\mathit{ij}}={\delta }_{i+N+1\text{,}j}\text{,}i=0\text{,}1\text{,}\dots \text{,}N(N+1)-1\text{,}I$ is an (N + 1) × (N + 1) identity matrix and δ denote the kronecker delta;

2. ${\mathit{tu}}_N(x\text{,}t)=\underline{u}{\hat{\mu }}_t{\underline{\varphi }}_{x\text{,}t}$ , where ${\hat{\mu }}_t=\Phi {\tilde{\mu }}_t{\Phi }^{-1}$ and ${\tilde{\mu }}_t=I\otimes \mu$ with the elements ${({\tilde{\mu }}_t)}_{\mathit{ii}}=\mu \text{,}i=0\text{,}1\text{,}\dots \text{,}N$ ;

3. $\frac{{\partial }^p}{\partial x^p}{u}_N(x\text{,}t)=\underline{u}{\hat{\eta }}_x^p{\underline{\varphi }}_{x\text{,}t}$ for $p\in \mathbb{N}$ , where ${\hat{\eta }}_x=\Phi {\tilde{\eta }}_x{\Phi }^{-1}$ and ${\tilde{\eta }}_x=\eta \otimes I$ with the elements ${({\tilde{\eta }}_x)}_{\mathit{ij}}=(j+1){\delta }_{i\text{,}j+1}I\text{,}j=0\text{,}1\text{,}\dots \text{,}N-1$ ;

4. $\frac{{\partial }^q}{\partial t^q}{u}_N(x\text{,}t)=\underline{u}{\hat{\eta }}_t^q{\underline{\varphi }}_{x\text{,}t}$ for $q\in \mathbb{N}$ , where ${\hat{\eta }}_t=\Phi {\tilde{\eta }}_t{\Phi }^{-1}$ and ${\tilde{\eta }}_t=I\otimes \eta$ with the elements ${({\tilde{\eta }}_t)}_{\mathit{ii}}=\eta \text{,}i=0\text{,}1\text{,}\dots \text{,}N$ .

1. By Corollary 2.2 and kronecker product properties, we can write $$\mathit{xu}(x\text{,}t)=x(\underline{u}{\underline{\varphi }}_{x\text{,}t})=\underline{u}\Phi x({\underline{X}}_x\otimes {\underline{X}}_t)=\underline{u}\Phi (x{\underline{X}}_x)\otimes {\underline{X}}_t=\underline{u}\Phi {\tilde{\mu }}_x{\Phi }^{-1}{\underline{\varphi }}_{x\text{,}t}=\underline{u}{\hat{\mu }}_x{\underline{\varphi }}_{x\text{,}t}\text{,}$$ where ${\tilde{\mu }}_x=\mu \otimes I$ and so $${({\tilde{\mu }}_x)}_{\mathit{ij}}=\left\{\begin{array}{cc}1\text{,}\hfill & j=i+N+1\text{,}i=0\text{,}1\text{,}\dots \text{,}N(N+1)-1\text{,}\hfill \\ 0\text{,}\hfill & \text{otherwise.}\hfill \end{array}\right.$$ which can be written as ${({\tilde{\mu }}_x)}_{\mathit{ij}}={\delta }_{i+N+1\text{,}j}\text{,}i=0\text{,}1\text{,}\dots \text{,}N(N+1)-1$ .

2. The proof is similar to the proof of part (a).

3. By Corollary 2.2 and kronecker product properties, we have $$\frac{{\partial }^p}{\partial x^p}{u}_N(x\text{,}t)=\frac{{\partial }^p}{\partial x^p}(\underline{u}{\underline{\varphi }}_{x\text{,}t})=\underline{u}\Phi \frac{{\partial }^p}{\partial x^p}({\underline{X}}_x\otimes {\underline{X}}_t)=\underline{u}\Phi \left(\frac{{\partial }^p}{\partial x^p}{\underline{X}}_x\right)\otimes {\underline{X}}_t=\underline{u}\Phi {\tilde{\eta }}_x^p{\Phi }^{-1}{\underline{\varphi }}_{x\text{,}t}=\underline{u}{\hat{\eta }}_x^p{\underline{\varphi }}_{x\text{,}t}\text{,}$$ where ${({\tilde{\eta }}_x)}_{\mathit{ij}}=(j+1){\delta }_{i\text{,}j+1}I\text{,}j=0\text{,}1\text{,}\dots \text{,}N-1$ .

4. The proof is similar to the proof of part (c). □

The effect of $x^p{t}^q\text{,}p\text{,}q\in \mathbb{N}$ on the coefficients of u_N(x,t) = u ϕ_x,t is equivalent to post multiplication of u by $\widehat{{\mu }_x^p{\mu }_t^q}$ , i.e. $$x^p{t}^q{u}_N(x\text{,}t)=\underline{u}\widehat{{\mu }_x^p{\mu }_t^q}{\underline{\varphi }}_{x\text{,}t}\mathit{or}{x}^p{t}^q{\underline{\varphi }}_{x\text{,}t}=\widehat{{\mu }_x^p{\mu }_t^q}{\underline{\varphi }}_{x\text{,}t}\text{,}$$ where $\widehat{{\mu }_x^p{\mu }_t^q}=\Phi {\tilde{\mu }}_x^p{\tilde{\mu }}_t^q{\Phi }^{-1}$ and the elements of the matrices ${\tilde{\mu }}_x^p$ and ${\tilde{\mu }}_t^q$ are determined by $${\left({\tilde{\mu }}_x^p\right)}_{\mathit{ij}}=\left\{\begin{array}{cc}1\text{,}\hfill & j=i+p(N+1)\text{,}i=0\text{,}1\text{,}\dots \text{,}(N-p+1)(N+1)-1\text{,}\hfill \\ 0\text{,}\hfill & \text{otherwise}\hfill \end{array}\right.$$ and $${\left({\tilde{\mu }}_t^q\right)}_{\mathit{ij}}=\left\{\begin{array}{cc}{\mu }^q\text{,}\hfill & i=j\text{,}i=0\text{,}1\text{,}\dots \text{,}N\text{,}\hfill \\ 0\text{,}\hfill & \text{otherwise}\hfill \end{array}\right.$$ with (μ^q)_ij = δ_i+q,j,i = 0,1, … ,N − q.

In addition, the matrix ${\tilde{\mu }}_x^p{\tilde{\mu }}_t^q$ has the following simple form $${\tilde{\mu }}_x^p{\tilde{\mu }}_t^q={\tilde{\mu }}_t^q{\tilde{\mu }}_x^p=\left(\begin{array}{cc}\hfill \overline{0}\hfill & \hfill B\hfill \\ \hfill \overline{0}\hfill & \hfill \overline{0}\hfill \end{array}\right)\text{,}$$ where $\overline{0}$ is an p × p zero matrix and B is an (N − p + 1) × (N − p + 1) diagonal matrix with the elements μ^q.

The matrices ${\tilde{\eta }}_x^p$ and ${\tilde{\eta }}_t^q$ have the following simple forms $${\tilde{\eta }}_x^p=\left(\begin{array}{cc}\hfill \overline{0}\hfill & \hfill \overline{0}\hfill \\ \hfill C\hfill & \hfill \overline{0}\hfill \end{array}\right)\text{,}$$ where $$C=\mathit{diag}((p(p-1)\times \cdots \times 1)I\text{,}((p+1)p\times \cdots \times 2)I\text{,}\dots \text{,}(N(N-1)\times \cdots \times (N-P+1))I)\text{,}$$ and $${\tilde{\eta }}_t^q=\mathit{diag}({\eta }^q\text{,}{\eta }^q\text{,}\dots \text{,}{\eta }^q)\text{.}$$

Let u_N(x,t) = u ϕ_x,t =  uΦX_x,t and v_N(x,t) = v ϕ_x,t =  vΦX_x,t, where v is a vector similar to u with elements v_i. Then $$u_N(x\text{,}t)v_N(x\text{,}t)=\underline{u}\widehat{V}{\underline{\varphi }}_{x\text{,}t}=\underline{u}\Phi V{\underline{X}}_{x\text{,}t}\text{,}$$ where $\widehat{V}=\Phi V{\Phi }^{-1}$ and Φ_j’s are the columns of matrix Φ with $$V=\left(\begin{array}{cccc}\hfill \underline{v}{\Phi }_0\hfill & \hfill \underline{v}{\Phi }_1\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_{{(N+1)}^2-1}\hfill \\ \hfill 0\hfill & \hfill \underline{v}{\Phi }_0\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_{{(N+1)}^2-2}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_0\hfill \end{array}\right)\text{.}$$

By assumptions of lemma, we have $$u_N(x\text{,}t)v_N(x\text{,}t)=\underline{u}\Phi ({\underline{X}}_{x\text{,}t}\times (\underline{v}\Phi {\underline{X}}_{x\text{,}t}))\text{,}$$ thus it suffices to show X_x,t × vΦX_x,t) = VX_x,t. By a simple computation, we can write $${\underline{X}}_{x\text{,}t}\underline{v}\Phi {\underline{X}}_{x\text{,}t}={[(1\text{,}t\text{,}\dots \text{,}{t}^N)\text{,}(1\text{,}t\text{,}\dots \text{,}{t}^N)x\text{,}\dots \text{,}(1\text{,}t\text{,}\dots \text{,}{t}^N)x^N]}^T(v_0\text{,}{v}_1\text{,}\dots \text{,}{v}_{{(N+1)}^2-1})({\Phi }_0|{\Phi }_1|\dots |{\Phi }_{{(N+1)}^2-1}){\underline{X}}_{x\text{,}t}={\left[(1\text{,}t\text{,}\dots \text{,}{t}^N)\underline{v}\text{,}(1\text{,}t\text{,}\dots \text{,}{t}^N)x\underline{v}\text{,}\dots \text{,}(1\text{,}t\text{,}\dots \text{,}{t}^N)x^N\underline{v}\right]}^T\times ({\Phi }_0|{\Phi }_1|\dots |{\Phi }_{{(N+1)}^2-1}){\underline{X}}_{x\text{,}t}={\left((\underline{v}{\Phi }_0){\underline{X}}_{x\text{,}t}|(\underline{v}{\Phi }_1){\underline{X}}_{x\text{,}t}|\dots |(\underline{v}{\Phi }_{{(N+1)}^2-1}){\underline{X}}_{x\text{,}t}\right)}^T{\underline{X}}_{x\text{,}t}=\left(\begin{array}{ccccc}\hfill \underline{v}{\Phi }_0\hfill & \hfill \underline{v}{\Phi }_1\hfill & \hfill \underline{v}{\Phi }_2\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_{{(N+1)}^2-1}\hfill \\ \hfill 0\hfill & \hfill \underline{v}{\Phi }_0\hfill & \hfill \underline{v}{\Phi }_1\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_{{(N+1)}^2-2}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \underline{v}{\Phi }_0\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_{{(N+1)}^2-3}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill \underline{v}{\Phi }_0\hfill \end{array}\right){\underline{X}}_{x\text{,}t}.\square$$

If u_N(x,t) = u ϕ_x,t, then $$u_N^k(x\text{,}t)=\underline{u}{\widehat{U}}^{k-1}{\underline{\varphi }}_{x\text{,}t}\text{,}k\in \mathbb{N}\text{,}$$ where $\widehat{U}=\Phi U^{k-1}{\Phi }^{-1}$ and U is an upper triangular matrix with elements $$U_{\mathit{ij}}=\left\{\begin{array}{cc}{\displaystyle \sum _{r=0}^{{(N+1)}^2-1}}{u}_r{\varphi }_{r\text{,}j-i}\text{,}\hfill & i⩽j\text{,}i=0\text{,}1\text{,}\dots \text{,}{(N+1)}^2-1\text{,}\hfill \\ 0\text{,}\hfill & i>j\text{.}\hfill \end{array}\right.$$

If u_N(x,t) = u ϕ_x,t, then

1. ${\left(\frac{{\partial }^p}{\partial x^p}{u}_N(x\text{,}t)\right)}^r=\underline{u}{\Psi }^p{M}_p^{r-1}{\Phi }^{-1}{\underline{\varphi }}_{x\text{,}t}$ for $p\text{,}r\in \mathbb{N}$ , where ${\Psi }^p=\Phi {\tilde{\eta }}_x^p$ and $$M_p=\left(\begin{array}{cccc}\hfill \underline{u}{\Psi }_0^p\hfill & \hfill \underline{u}{\Psi }_1^p\hfill & \hfill \cdots \hfill & \hfill \underline{u}{\Psi }_{{(N+1)}^2-1}^p\hfill \\ \hfill 0\hfill & \hfill \underline{u}{\Psi }_0^p\hfill & \hfill \cdots \hfill & \hfill \underline{u}{\Psi }_{{(N+1)}^2-2}^p\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill \underline{u}{\Psi }_0^p\hfill \end{array}\right)\text{.}$$

2. ${\left(\frac{{\partial }^q}{\partial t^q}{u}_N(x\text{,}t)\right)}^s=\underline{u}{\mathrm{\Gamma }}^q{N}_q^{s-1}{\Phi }^{-1}{\underline{\varphi }}_{x\text{,}t}$ for $q\text{,}s\in \mathbb{N}$ , where ${\mathrm{\Gamma }}^q=\Phi {\tilde{\eta }}_t^q$ and $$N_q=\left(\begin{array}{cccc}\hfill \underline{u}{\mathrm{\Gamma }}_0^q\hfill & \hfill \underline{u}{\mathrm{\Gamma }}_1^q\hfill & \hfill \cdots \hfill & \hfill \underline{u}{\mathrm{\Gamma }}_{{(N+1)}^2-1}^q\hfill \\ \hfill 0\hfill & \hfill \underline{u}{\mathrm{\Gamma }}_0^q\hfill & \hfill \cdots \hfill & \hfill \underline{u}{\mathrm{\Gamma }}_{{(N+1)}^2-2}^q\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill \underline{u}{\mathrm{\Gamma }}_0^q\hfill \end{array}\right)\text{.}$$

If u_N(x,t) = u ϕ_x,t, then

1. $u_N^k(x\text{,}t)\frac{{\partial }^p}{\partial x^p}{u}_N(x\text{,}t)=\underline{u}{\widehat{U}}^{k-1}{\widehat{M}}_p{\underline{\varphi }}_{x\text{,}t}\text{,}k\text{,}p\in \mathbb{N}$ ;

2. $u_N^k(x\text{,}t)\frac{{\partial }^q}{\partial t^q}{u}_N(x\text{,}t)=\underline{u}{\widehat{U}}^{k-1}{\widehat{N}}_q{\underline{\varphi }}_{x\text{,}t}\text{,}k\text{,}q\in \mathbb{N}$ ,

where ${\widehat{U}}^{k-1}=\Phi U^{k-1}{\Phi }^{-1}\text{,}{\widehat{M}}_p=\Phi M_p{\Phi }^{-1}$ and ${\widehat{N}}_q=\Phi N_q{\Phi }^{-1}$ .

The matrix ${\widehat{\Pi }}_D$ in Eq. (13) has the following structure $${\widehat{\Pi }}_D=\rho c_p{\hat{\eta }}_t-k_0\left[{\hat{\eta }}_x+{\alpha }_1(\Psi M_1{\Phi }^{-1}+{\widehat{M}}_2)+{\alpha }_2(2\Phi {\overline{M}}_1{\Phi }^{-1}+\widehat{U}{\widehat{M}}_2)\right]\text{,}$$ where $${\overline{M}}_1=\left(\begin{array}{cccc}\hfill \underline{u}{({\stackrel{\sim }{M}}_1)}_0\hfill & \hfill \underline{u}{({\stackrel{\sim }{M}}_1)}_1\hfill & \hfill \cdots \hfill & \hfill \underline{u}{({\stackrel{\sim }{M}}_1)}_{{(N+1)}^2-1}\hfill \\ \hfill 0\hfill & \hfill \underline{u}{({\stackrel{\sim }{M}}_1)}_0\hfill & \hfill \cdots \hfill & \hfill \underline{u}{({\stackrel{\sim }{M}}_1)}_{{(N+1)}^2-2}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill \underline{u}{({\stackrel{\sim }{M}}_1)}_0\hfill \end{array}\right)\text{,}$$ and ${\stackrel{\sim }{M}}_1=\Psi M_1$ .

Consider the following nonlinear transient heat conduction equation $$\frac{\partial }{\partial x}\left(k(u)\frac{\partial u}{\partial x}\right)+f(x\text{,}t)=\frac{\partial u}{\partial t}\text{,}x\in [0\text{,}1]\text{,}t>0\text{,}$$ where $$k(u)=1+u^2\text{,}$$ and $$f(x\text{,}t)=1-2x-2t\text{,}$$ with the supplementary conditions $$u(0\text{,}t)=t\text{,}u(1\text{,}t)=1+t\text{,}u(x\text{,}0)=x\text{.}$$ The exact solution of this problem is u(x,t) = x + t. For N = 2 and using the standard basis, the proposed method gives the nonlinear system $$\left\{\begin{array}{c}{u}_0=u_2=u_6=0\text{,}\hfill \\ u_1=u_3=1\text{,}\hfill \\ u_4+u_7=0\text{,}\hfill \\ u_5+u_8=0\text{,}\hfill \\ u_0(2u_0{u}_7+2u_1{u}_6+4u_3{u}_4)+u_1(2u_0{u}_6+2u_3^2)-2u_2+2u_7=2\text{,}\hfill \\ u_0(2u_0{u}_8+2u_1{u}_7+2u_2{u}_6+4u_3{u}_5+2u_4^2)\hfill \\ +u_1(2u_0{u}_7+2u_1{u}_6+4u_3{u}_4)\hfill \\ +u_2(2u_0{u}_6+2u_3^2)+2u_8=0\text{,}\hfill \end{array}\right.$$ with the solution {u₁ = u₃ = 1, u_i = 0 for i ≠ 1,3} which leads to the exact solution. Indeed, as mentioned in Ortiz and Samara (1981), the operational Tau method for equations with polynomial solution is exact, whenever the degree of the Tau approximation is at least equal to the degree of solution.

Consider a plane wall having variable thermo-physical properties (thermal conductivity, specific heat and density) in which its surface temperatures are related to the time. We assume that a heat sink or heat source is presented in the wall in which the magnitude of released or dissipated energy is a nonlinear function of space and temperature (for examples chemical reaction or electrical resistance). Also, the ambient temperature is zero. The following nonlinear heat conduction equation can be obtained using the first law of thermodynamics $$\frac{\partial }{\partial x}\left(k(u)\frac{\partial u}{\partial x}\right)+f(x\text{,}t)=u^2\frac{\partial u}{\partial t}\text{,}x\in [0\text{,}1]\text{,}t>0\text{.}$$ In a special case we assume that thermo-physical properties and sum of heat flux are $$k(u)=1+u+\frac{1}{2}{u}^2\text{,}$$ and $$f(x\text{,}t)=-{\sin }^{2}t(1+x\sin t-x^3\cos t)\text{.}$$ The supplementary conditions are assumed to be of the Dirichlet kind, namely $$u(0\text{,}t)=0\text{,}u(1\text{,}t)=\sin t\text{,}u(x\text{,}0)=0\text{.}$$ The exact solution of this problem is u(x, t) = x sin t.

Table 1 shows the absolute errors with respect to shifted Chebyshev basis functions at the selected points for various choices of N.

Consider the nonlinear transient heat conduction equation of the form $$\frac{\partial }{\partial x}\left(k(u)\frac{\partial u}{\partial x}\right)+f(x\text{,}t)=\frac{\partial u}{\partial t}\text{,}x\in [0\text{,}1]\text{,}t>0\text{.}$$ with supplementary conditions $$u(0\text{,}t)=0\text{,}u(1\text{,}t)+\frac{\partial u(1\text{,}t)}{\partial x}=2e^t\text{,}u(x\text{,}0)=x\text{,}$$ where thermo-physical properties and sum of heat flux are $$k(u)=1+u^2\text{,}$$ and $$f(x\text{,}t)={\mathit{xe}}^t(1-2e^{2t})\text{.}$$ The exact solution of the problem is u(x,t) = xe^t.

The absolute and estimate errors for this example with respect to shifted Chebyshev basis for various choices of N, reported in Tables 2 and 3, show the accuracy and efficiency of the proposed method. These results confirm that the absolute and estimate errors are in good agreement.

To check the stability of the proposed method, we perturb the coefficients of approximate solution by ε = 10⁻³, 10⁻⁵ and 10⁻⁷. Then, we solve the perturbed problem by the method and find out that there are no total changes in the final results. Table 4 shows the maximum absolute errors of the perturbed problem with respect to shifted Chebyshev basis for various choices of N.