Lanczos Tau method

In mathematics, the Lanczos tau method is a spectral method for the numerical solution of differential equations, named after Cornelius Lanczos, who first presented it in 1938. Unlike the Galerkin method, which requires the basis functions to satisfy the boundary conditions a priori, the tau method uses a set of orthogonal polynomials as basis functions. In the Tau method approach, the boundary conditions are imposed on the linear combination of basis functions, while the residual is made orthogonal to as many basis functions as possible — rather than to the entire set.

Overview

The tau method seeks an approximation to the solution of a differential equation by expanding the solution in terms of a set of basis functions, typically orthogonal polynomials such as Chebyshev or Legendre polynomials. The key idea is to require that the residual of the differential equation be orthogonal to as many basis functions as possible, while the boundary or auxiliary conditions are forced to be exactly satisfied.

Consider a differential equation of the form

${\displaystyle Dy(x)=f(x),\quad x\in (a,b)}$

subject to ${\textstyle \nu }$ appropriate boundary or auxiliary conditions

${\displaystyle g_{i}(y)=s_{i},\quad i=1,\ldots ,\nu .}$

Here, ${\textstyle D}$ is a differential operator that can be, at least approximately, written in the form

${\displaystyle Dy=\sum _{i=0}^{\nu }u_{i}{\dfrac {d}{dx^{i}}}y,}$

${\textstyle \nu }$ is the differential problem order, ${\displaystyle f,g_{i},s_{i},i=1,\ldots ,\nu }$ and ${\textstyle u_{i},i=0,\ldots ,\nu }$ are given. For the tau method to be applied, it is assumed that the solution admits a generalised Fourier series

${\displaystyle y(x)=\sum _{j=0}^{\infty }a_{j}P_{j}(x)}$

where ${\displaystyle {\mathcal {P}}=\{P_{j}\}_{j=0}^{\infty }}$ is a set of square-integrable orthogonal functions defined on ${\displaystyle [a,b]}$, and ${\displaystyle a_{j}}$ are unknown coefficients. Usually ${\displaystyle {\mathcal {P}}}$ is an orthogonal polynomial basis.

The Tau Formulation

In the tau method, for some integer ${\displaystyle n>\nu }$, an approximation to ${\displaystyle y}$ is found in the subspace ${\displaystyle \mathbb {P} _{n}}$ spanned by ${\displaystyle \{P_{j}\}_{j=0}^{n}}$, in the form

${\displaystyle y_{n}(x)=\sum _{j=0}^{n}a_{j}^{(n)}P_{j}(x)}$

where the coefficients ${\displaystyle a_{j}^{(n)}}$ are found by requiring that the boundary conditions hold

${\displaystyle g_{i}(y_{n})=s_{i},\quad i=1,\ldots ,\nu }$

and using the remaining ${\displaystyle n+1-\nu }$ degrees of freedom to impose that ${\displaystyle Dy_{n}}$ agrees with ${\displaystyle Dy}$ as far as possible.The last condition means that the residual is orthogonal to each ${\displaystyle P_{j},\ j=0,\ldots ,n-\nu }$ or, equivalently, forces the residual to be an order ${\displaystyle n+1-\nu }$ Fourier series in ${\displaystyle {\mathcal {P}}}$

${\displaystyle \tau \equiv Dy-Dy_{n}=\sum _{j=n+1-\nu }^{\infty }\tau _{j}P_{j}}$

In case all coefficients ${\displaystyle u_{i}}$ in ${\displaystyle D}$ are polynomials, and if ${\displaystyle d_{i}}$is the degree of ${\displaystyle u_{i}}$, then the residual ${\displaystyle \tau }$ is a polynomial of degree ${\displaystyle n+\max\{d_{i}-i,\ i=0,\ldots ,\nu \}}$.

Operational Approach of the Tau Method

In 1981, Ortiz and Samara proposed an operational approach to the Tau method based on operational matrices of differential operators. For a fixed ${\textstyle n>\nu }$, let ${\displaystyle {\mathcal {P}}_{n}=[P_{0},P_{1},\ldots ,P_{n}]}$ be a basis for the subspace ${\displaystyle \mathbb {P} _{n}}$, and let ${\displaystyle {\mathsf {a}}_{n}=[a_{0}^{(n)},a_{1}^{(n)},\ldots ,a_{n}^{(n)}]^{T}}$ be the vector of the Fourier coefficients of ${\displaystyle y_{n}}$, then the approximate solution can be written as ${\textstyle y_{n}={\mathcal {P}}_{n}{\mathsf {a}}_{n}}$.

Let ${\displaystyle x{\mathcal {P}}_{n}=[xP_{0},xP_{1},\ldots ,xP_{n}]}$ denote the multiplication of basis elements by ${\displaystyle x}$, and let ${\displaystyle {\mathcal {P}}'_{n}=[P'_{0},P'_{1},\ldots ,P'_{n}]}$ denote differentiation of basis elements. These operations can be represented as matrix operations by

${\displaystyle x{\mathcal {P}}_{n}={\mathcal {P}}_{n}{\mathsf {M}}\quad {\text{and}}\quad {\mathcal {P}}'_{n}={\mathcal {P}}_{n}{\mathsf {N}}}$

Matrices ${\textstyle {\mathsf {M}}}$ and ${\textstyle {\mathsf {N}}}$ can be evaluated from properties of ${\textstyle {\mathcal {P}}_{n}}$. If ${\textstyle {\mathcal {P}}_{n}}$ is an orthogonal polynomial sequence, then a three-term recurrence relation holds:

${\displaystyle {\begin{cases}xP_{j}=\alpha _{j}P_{j+1}+\beta _{j}P_{j}+\gamma _{j}P_{j-1},&j=0,1,\ldots ,n\\P_{0}=1;\quad P_{-1}=0\end{cases}}.}$

Then ${\textstyle {\mathsf {M}}}$ is the Jacobi matrix ${\textstyle {\mathsf {M}}=[\mu _{i,j}]_{i,j=0}^{n+1,n}}$ with

${\displaystyle \mu _{i,j}={\begin{cases}\gamma _{j},&i=j-1,\\\beta _{j},&i=j,\\\alpha _{j},&i=j+1,\\0,&|i-j|>1\end{cases}}}$

and ${\textstyle {\mathsf {N}}=[\eta _{i,j}]_{i,j=0}^{n,n}}$ where ${\textstyle \eta _{i,j}}$ are the connection coefficients in

${\displaystyle P'_{j}=\sum _{i=0}^{j-1}\eta _{i,j}P_{i}}$

From matrices ${\textstyle {\mathsf {M}}}$ and ${\textstyle {\mathsf {N}}}$, the action of the differential operator over polynomials is translated to the action of the algebraic operator, over the vector of Fourier coefficients, defined by the operational matrix

${\displaystyle {\mathsf {D}}=\sum _{i=0}^{\nu }u_{i}({\mathsf {M}}){\mathsf {N}}^{i},}$

in the sense that ${\displaystyle {\mathsf {D}}{\mathsf {a}}_{n}={\mathsf {f}}}$, where ${\displaystyle {\mathsf {f}}=[f_{0},f_{1},\ldots ,f_{n}]}$ is the vector of Fourier coefficients of the projection of ${\displaystyle f}$ onto ${\displaystyle \mathbb {P} _{n}}$. If, for the coefficients of ${\displaystyle D}$, we have

${\displaystyle u_{i}=\sum _{j=0}^{d_{i}}u_{i,j}P_{j}}$

at least approximately, then

${\displaystyle u_{i}({\mathsf {M}})=\sum _{j=0}^{d_{i}}u_{i,j}P_{j}({\mathsf {M}}).}$

The Tau matrix

The action of boundary conditions over the base ${\textstyle {\mathcal {P}}_{n}}$ is represented by the ${\displaystyle \nu \times (n+1)}$ matrix ${\displaystyle {\mathsf {G}}=[g_{i,j}]_{i,j=0}^{\nu -1,n}}$ with

${\displaystyle g_{i,j}=g_{i+1}(P_{j})}$

and the vector ${\displaystyle {\mathsf {s}}=[s_{1},s_{2},\ldots ,s_{i}]^{T}}$. Then ${\textstyle y_{n}={\mathcal {P}}_{n}{\mathsf {a}}_{n}}$, is an approximate solution to

${\displaystyle {\begin{cases}Dy(x)=f(x),&x\in (a,b)\\g_{i}(y)=s_{i},&i=1,\ldots ,\nu \end{cases}}}$

in the Tau sense if

${\displaystyle {\mathsf {T}}_{n}{\mathsf {a}}_{n}={\mathsf {b}}_{n}}$

where

${\displaystyle {\mathsf {T}}_{n}=\left[{\begin{array}{l}{\mathsf {G}}\\{\mathsf {D}}_{n}\end{array}}\right],\quad {\text{and}}\quad {\mathsf {b}}_{n}=\left[{\begin{array}{l}{\mathsf {s}}\\{\mathsf {f}}_{n}\end{array}}\right],}$

${\displaystyle {\mathsf {D}}_{n}}$ is matrix ${\displaystyle {\mathsf {D}}}$ truncated to its first ${\displaystyle n+1-\nu }$ rows and ${\displaystyle {\mathsf {f}}_{n}}$ is vector ${\displaystyle {\mathsf {f}}}$ truncated to its first ${\displaystyle n+1-\nu }$ elements.

See also

Galerkin method

References