---
title: 11-6. Spectral integration and two point BVP
tags: Final reports
---
**Reference**
1. Spectral Integration and Two-Point Boundary Value Problems, L. Greengard, SIAM Journal of Numerical Analysis **28** (1991) 1071–1080. [10.1137/0728057](https://doi.org/10.1137/0728057)
---
### Spetral integration and two-point boundary value problems(亞衛)
#### 1. Chebyshev approxamation by finding coefficients
This article uses the first kind of Chebyshev polynomial.
The Chebyshev polynomial of degree k is defined by
$$ T_k(\cos\theta) = \cos(k\theta) $$ with recursive relation$$T_0(x) = 1,\quad T_1(x)=x, \\
T_{k+1}(x)=2xT_k(x)-T_{k-1}(x) \quad\text{for }k\geq1.$$
This article also uses the first kind of Chebyshev nodes defined as
$x_j = \cos\left( \dfrac{(2j+1)\pi}{2N}\right),\quad j=0,1,...,N-1.$
Therefore write $u(x)$ in Chebyshev series, we have
\begin{aligned}
u(x)
&=\dfrac{1}{2}a_0T_0(x)+a_1T_1(x)+a_2T_2(x)+\cdots \\ &={\sum_{k=0}^{\infty}}{^\prime}{a_kT_k(x)}
\end{aligned}
with coefficients given by
$$a_k=\dfrac{2}{\pi}\int^1_{-1}f(x)T_k(x)(1-x^2)^{-1/2}dx$$ or
$$a_k = \dfrac{2}{N}\sum_{j=0}^{N-1} f(x_j)T_k(x_j).$$
* __Differentiation__
Write $u'(x)$ in the Chebyshev series expansion with coefficients $\{b_k\}$:
$$u^{\prime}(x)={\sum_{k=0}^{N}}{^\prime}{b_kT_k(x)}$$
Then the coefficients $\{b_k\}$ has the recursive relation with $\{a_k\}$ as follow:
$$b_k=\sum_{p\ =\ k+1\\ p+k \ \ odd}^{\infty} p\cdot a_p.$$
* __Integration__
Write $\int_{-1}^{x} u(t)dt$ in the Chebyshev series expansion with coefficients $\{d_k\}$
$$\int^x_{-1}u(t)dt=
{\sum_{k=0}^{N}}{^\prime}{d_kT_k(x)}$$
Then the coefficients $\{d_k\}$ has the recursive relation with $\{a_k\}$ as follow:
$$d_k=\dfrac{1}{2k}(a_{k-1}-a_{k+1}) \quad\text{for}\quad k\geq 1\\
d_0=d_1-d_2+d_3-d_4+\cdots .$$
* __Error__
The process of differentiation via Chebyshev series amplifies error by a factor proportional to $N^2$.
Now let
$$\hat a = (a_1+\epsilon,a_2+\epsilon,\cdots),\quad \hat b = \mathscr{D} (\hat a).$$
where $\epsilon$ is a small number, then
$$ \dfrac{||\hat b-b||_{\infty}}{||\hat a-a||_{\infty}} = O(N^2).$$
But the process of integration via Chebyshev series amplifies error by a factor of less than 2.4, i.e.,
$$ \dfrac{|\hat{d_k}-d_k|} {||\hat a-a||_{\infty}} \leq 1\quad\text{ for } k\geq1$$
and
$$ \dfrac{|\hat{d_0}-d_0|} {||\hat a-a||_{\infty}} \leq \frac{1}{2}+\frac{1}{4}+\zeta(2) < 2.4$$
where $\zeta$ denotes the Reimann zeta function and $\hat d = \mathscr{I} (\hat a)$.
---
#### 2. Differentiation matrix and Integration matrix via Chebyshev polynomial
Since the chebyshev polynomial of degree $k$ is defined as $T_k(\cos\theta) = \cos(k\theta)$, using the Chebtshev points of the fisrt kind $x_j = \cos\left( \dfrac{(2j+1)\pi}{2N}\right)$, we have $$ T_k(x_j)= T_k(\cos\left( \dfrac{(2j+1)\pi}{2N}\right)) = \cos\left( \dfrac{(2j+1)\cdot k \pi}{2N}\right).$$
Let $$ \theta_j = \dfrac{(2j+1)\pi}{2N},\quad j=0,1,\cdots,N-1 $$ gives
$$T_k(x_j) =T_k(\cos(\theta_j))= \cos(k\cdot \theta_j),\quad\text{for } j=0,1,\cdots,N-1.$$
Now since
\begin{aligned}
a_k & = \dfrac{2}{N}\sum_{j=0}^{N-1} f(x_j)T_k(x_j) \\
&=\dfrac{2}{N}\sum_{j=0}^{N-1} f(x_j)\cos(k\cdot\theta_j)
\end{aligned}
Write $\{f(x_j)\}$ and $\{a_k\}$ into matrices form:
$$\dfrac{2}{N}
\begin{pmatrix}
\cos(0) & \cos(0) & \cdots & \cos(0) \\
\cos(\theta_0) & \cos(\theta_1) & \cdots & \cos(\theta_{N-1}) \\
\cos(2\theta_0)& \cos(2\theta_1) & \cdots & \cos(2\theta_{N-1})\\
\vdots & \vdots& \ddots & \vdots \\
\cos((N-1)\theta_0) & \cos((N-1)\theta_1) & \cdots &\cos((N-1)\theta_{N-1})\\
\end{pmatrix}
\begin{pmatrix}
f(x_0) \\ f(x_1) \\f(x_2)\\ \vdots \\ f(x_{N-1})
\end{pmatrix}=
\begin{pmatrix}
a_0 \\ a_1 \\a_2\\ \vdots \\ a_{N-1}
\end{pmatrix}
$$
or
$$ CF={\bf a}
$$
where $C_{N\times N}$ is the transformation matrix with $c_{ij} = \frac{2}{N}\cos((i-1)\theta_{j-1}).$
Moreover, if we want to transfrom the coefficients $\{a_k\}$ back to the function value $\{f_k\}$, we just have to times a inverse function of transformation matrix $C^{-1}$, i.e.,
$$F=C^{-1}{\bf a}.$$
* __Differentiation matrix__
Since the coefficients $\{b_k\}$ has the recursive relation with $\{a_k\}$ as follows:
$$b_k=\sum_{p\ =\ k+1\\ p+k \ \ odd}^{\infty} p\cdot a_p.$$
More percisely,
\begin{aligned}
b_0 &= 1\cdot a_1 + 3\cdot a_3 + 5\cdot a_5 + \cdots\\
b_1 &= 2\cdot a_2 + 4\cdot a_4 + 6\cdot a_6+\cdots\\
b_2 &= 3\cdot a_3 + 5\cdot a_5 + 7\cdot a_7+\cdots\\
&\quad\vdots
\end{aligned}
write $\{b_k\}$ and $\{a_k\}$ into matrices form:
$$
\begin{pmatrix}
0 & 1 & 0 & 3 & 0 & 5 & \cdots \\
0 & 0 & 2 & 0 & 4 & 0 & \cdots \\
0 & 0 & 0 & 3 & 0 & 5 & \cdots \\
\vdots & \vdots& &\ddots & \ddots & \vdots & \vdots\\
\vdots & \vdots &&& \ddots &&\vdots\\
\vdots & \vdots &&&& \ddots &\vdots\\
0 & 0 & \cdots& \cdots &\cdots & \cdots &0\\
\end{pmatrix}
\begin{pmatrix}
a_0 \\ a_1 \\a_2\\ \vdots \\ \vdots \\ \vdots \\ a_{N-1}
\end{pmatrix} =
\begin{pmatrix}
b_0 \\ b_1 \\b_2\\ \vdots \\ \vdots \\ \vdots \\b_{N-1}
\end{pmatrix}
$$
or
$$D{\bf a}={\bf b}$$
where $D_{N\times N}$ is the __differentiation matrix__ which transforms coefficients $\{a_k\}$ to $\{b_k\}$.
Note that the first column and the last row of $D$ are all zeros.
However we can rewrite matrix $D$ in a smaller size:
$$
\begin{pmatrix}
1 & 0 & 3 & 0 & 5 & \cdots \\
0 & 2 & 0 & 4 & 0 & \cdots \\
0 & 0 & 3 & 0 & 5 & \cdots \\
\vdots & \vdots & & \ddots && \vdots \\
\vdots & \vdots & && \ddots & \vdots \\
\end{pmatrix}
\begin{pmatrix}
a_1 \\ a_2 \\a_3\\ \vdots \\ \vdots \\ a_{N-1}
\end{pmatrix} =
\begin{pmatrix}
b_0 \\ b_1 \\b_2\\ \vdots \\ \vdots \\b_{N-2}
\end{pmatrix}
$$
or
$$D'{\bf a}={\bf b}$$
where $D'_{(N-1) \times (N-1)}$ is the deduced differentiation matrix.
* __Integration matrix__
Since the coefficients $\{d_k\}$ has the recursive relation with $\{a_k\}$ as follows:
$$d_k=\dfrac{1}{2k}(a_{k-1}-a_{k+1}) \quad\text{for}\quad k\geq 1\\
d_0=d_1-d_2+d_3-d_4+\cdots .$$
More percisely,
\begin{aligned}
d_1 &= \frac{1}{2}(a_0-a_2)=\frac{1}{2}a_0-\frac{1}{2}a_2\\
d_2 &= \frac{1}{4}(a_1-a_3)=\frac{1}{4}a_1-\frac{1}{4}a_3\\
d_3 &= \frac{1}{6}(a_2-a_4)=\frac{1}{6}a_2-\frac{1}{6}a_4\\
&\quad\vdots
\end{aligned}
Specially,
\begin{aligned}
d_0 &= d_1-d_2+d_3-d_4+d_5-d_6+\cdots\\
&= (\frac{1}{2}a_0-\frac{1}{2}a_2)+(\frac{1}{4}a_1-\frac{1}{4}a_3)+(\frac{1}{6}a_2-\frac{1}{6}a_4)+\cdots\\
&= \frac{1}{2}a_0-\frac{1}{4}a_1-\frac{1}{3}a_2+\frac{1}{8}a_3-\frac{1}{15}a_4+\cdots\\
&= \frac{1}{2}a_0-\frac{1}{4}a_1
-\frac{1}{1\times3}a_2
+\frac{1}{2\times4}a_3
-\frac{1}{3\times5}a_4+\cdots \\
&=\frac{1}{2}a_0-\frac{1}{4}a_1+\sum_{k=2}^{\infty}\dfrac{(-1)^{k+1}}{(k-1)(k+1)}\cdot a_k.
\end{aligned}
Write $\{d_k\}$ and $\{a_k\}$ into matrices form:
$$
\begin{pmatrix}
\dfrac{1}{2} & \dfrac{-1}{4} & \dfrac{-1}{1\times3} & \dfrac{1}{2\times4} & \dfrac{-1}{3\times5} &\cdots \\
\dfrac{1}{2} & 0 &\dfrac{-1}{2} \\
&\dfrac{1}{4} & 0 &\dfrac{-1}{4} \\
&&\dfrac{1}{6} & 0 &\dfrac{-1}{6} \\
&&&&\ddots
\end{pmatrix}
\begin{pmatrix}
a_0 \\ a_1 \\a_2\\ \vdots \\ \vdots \\ \vdots \\ a_{N-1}
\end{pmatrix} =
\begin{pmatrix}
d_0 \\ d_1 \\d_2\\ \vdots \\ \vdots \\ \vdots \\d_{N-1}
\end{pmatrix}
$$
or
$$I{\bf a}={\bf d}$$
where $I_{N\times N}$ is the __integration matrix__ which transforms coefficients $\{a_k\}$ to $\{d_k\}$.
In particular, the __spectral differentiation matrix __$D_N$ for Chebyshev nodes(od the fisrt kind) can be expresses in terms of $D$ by the formula
$$D_N=C^{-1}\cdot D\cdot C$$
where $C$ is the transformation matrix mentioned above;
On the other hand, the __spectral integration matrix__ $I_N$ for Chebyshev nodes(of the first kind) is defined by the formula
$$I_N=C^{-1}\cdot D\cdot C.$$
First try $u(x)=\sin(x)$:
```C++=
#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<cmath>
#include<vector>
#include<conio.h>
#include<complex>
#include<valarray>
using namespace std;
#define N 16
#define pi 3.14159265358979323846
double infinity_norm(const vector<double> &a) {
double max=0;
double k;
int n = a.size();
for(int i = 0; i < n; i++)
{
k = fabs(a[i]);
if(k > max) max = k;
}
return max;
}
void Peek(vector<double> vector) { //Output the vector
for (int i = 0; i < vector.size(); i++) {
cout <<"vec["<< i << "] = " << vector[i] << endl;
}
}
vector<double> Chebval(vector<double> coef,vector<double> x,int type){
//type 1: regular sum, no need scalar 0.5 on the first term
//type 2: Sigma prime, need scalar 0.5 on the first term
int k = coef.size();
vector<double> V(k,0.0);
for(int i=0;i<k;i++){
V[i] = 0.0;
for(int j=0;j<k;j++){
if(j==0 && type==2) V[i] += 0.5*coef[j]*cos(j*acos(x[i]));
else{
V[i] += coef[j]*cos(j*acos(x[i]));
}
}
}
return V;
}
double Err(vector<double> num, vector<double> exa){
int k = num.size();
vector<double> error(k,0.0);
for(int i=0;i<k;i++){
error[i] = abs(num[i] - exa[i]);
}
return infinity_norm(error);
}
int main(void) {
cout.precision(12);
vector<double> xx(N+1,0.0);
vector<double> f(N+1,0.0);
for(int i=0;i<N;i++){
xx[i] = cos((2*i + 1.0)*pi / (2*N));
f [i] = sin(xx[i]); //f(x)=sin(x)+0.01 sin(10x)
}
vector<double> coef(N+1,0.0);
for(int j=0;j<N;j++){
coef[j] = 0.0;
for(int i=0;i<N;i++){
coef[j] = coef[j] + f[i] * cos(pi*j*(2*i+1)/(2*N));
}
}
for(int i=1;i<N;i++) coef[i] = 2*coef[i] / N;
vector<double> val(N+1,0.0);
vector<double> error(N+1,0.0);
vector<double> a(N+1,0.0),ahat(N+1,0.0);
vector<double> b(N+1,0.0),bhat(N+1,0.0);
vector<double> d(N+1,0.0);
double e = 0.00001;
// vector a is the chebyshev coefficients
for(int i=0;i<N;i++){
a[i] = coef[i];
ahat[i] = a[i] + e;
}
// vector b is the differentiation coefficients
for(int i=0;i<N;i++){
for(int p=(i+1);p<N;p++){
if((p+i)%2 == 1){
b[i] += p*a[p];
bhat[i] += p*ahat[p];
}
}
}
vector<double> e1(N+1,0.0),e2(N+1,0.0);
for(int i=0;i<N;i++){
e1[i] = abs(ahat[i] - a[i]);
e2[i] = abs(bhat[i] - b[i]);
}
double conv = infinity_norm(e2)/infinity_norm(e1);
cout<<"conv error = "<<conv;
// Calculate the error from numerical solution and exact solution
val = Chebval(coef,xx,2);
cout<<"The first error is "<<Err(val,f)<<endl;
// Calculate the differentiation error
vector<double> diff(N+1,0.0);
for(int i=0;i<=N;i++) {
diff[i] = cos(xx[i]);
b[i] *= 2.0;
}
val = Chebval(b,xx,2);
cout<<endl<<"The second error is "<<Err(val,diff)<<endl;
// Calculate the integral coeffiecient
d[0] = 0.0;
for(int i=N-1;i>0;i--){
d[i] = (a[i-1]-a[i+1]) / (2*i);
d[0] += pow(-1.0,i+1)*d[i];
}
vector<double> Inte(N+1,0.0);
for(int i=0;i<=N;i++) Inte[i] = cos(1) - cos(xx[i]);
val = Chebval(d,xx,1);
cout<<"The third error is "<<Err(val,Inte)<<endl;
}
```
$$ u(x) = \sin(x),\quad
\hat u(x) = \sum_{k=0}^{\infty}{^\prime}{a_kT_k(x)} ,\quad Err1=||u(x)-\hat u(x)||_{\infty} \\
u'(x) = \cos(x),\quad
\hat v(x) = \sum_{k=0}^{\infty}{^\prime}{b_kT_k(x)} ,\quad Err2=||u'(x)-\hat v(x)||_{\infty} \\
\int_{-1}^{x}u(t)dt = \cos(-1)-\cos(x),\quad
\hat f(x) = \sum_{k=0}^{\infty}{^\prime}{d_kT_k(x)} ,\quad \\
Err3=||\int_{-1}^{x}u(t)dt-\hat f(x)||_{\infty} ,\\
L=\dfrac{||\hat b-b||_{\infty}}{||\hat a-a||_{\infty}}$$
|$N$| $Err1$ | $Err2$ |$Err3$ | $L$| $L/N^2$|
|---|:------ |:------ |:----- |:----| :------ |
|$16$ |$2.775557e-15$|$2.637889e-13$|$1.830133e-15$|$64$|$0.25$|
|$32$ |$4.107825e-15$|$4.257705e-13$|$2.220446e-16$|$256$|$0.25$|$0.25$|
|$64$ |$5.884182e-15$|$1.115663e-11$|$2.012279e-15$|$1024$|$0.25$|
|$128$|$8.770761e-15$|$8.168676e-11$|$1.389947e-15$|$4096$|$0.25$|
|$512$|$6.972201e-14$|$5.455472e-9$|$3.776121e-14$|$65536$|$0.25$|
|$1024$|$1.624950e-13$|$1.982892e-8$|$4.281724e-15$|$262144$|$0.25$|
|$4096$|$6.961093e-13$|$2.198919e-6$|$1.245135e-13$|$4194304$|$0.25$|
Then try $u(x) = \sin(x)+0.01\sin(10x)$.
The code is same as above except:
```C++=
f [i] = sin(xx[i]) + 0.01*sin(10*xx[i]); //f(x)=sin(x)+0.01 sin(10x)
diff[i] = cos(xx[i]) + 0.1*cos(10*xx[i]);
Inte[i] = cos(1) + 0.001*cos(10) - cos(xx[i]) - 0.001*cos(10*xx[i]);
```
$$u(x) = \sin(x)+0.01\sin(10x), u'(x) = \cos(x)+0.1\cos(10x)\\
\int_{-1}^{x}u(t)dt = \cos(1)+0.001\cos(10) - \cos(x) - 0.001\cos(10x)\\
Err1 = ||u(x)-\sum_{k=0}^{\infty}{^\prime}{a_kT_k(x)}||_{\infty} \\
Err2 = ||u'(x)-\sum_{k=0}^{\infty}{^\prime}{b_kT_k(x)}||_{\infty} \\
Err3 = ||\int_{-1}^{x}u(t)dt-\sum_{k=0}^{\infty}{^\prime}{d_kT_k(x)}||_{\infty} \\
$$
|$N$| $Err1$ | $Err2$ |$Err3$ | $L$| $L/N^2$|
|---|:------ |:------ |:----- |:----| :-------|
|$16$ |$2.109423e-15$|$0.00303208274$|$3.863301e-6$|$64$|$0.25$|
|$32$ |$4.996003e-15$|$3.186340e-13$|$2.641984e-15$|$256$|$0.25$|
|$64$ |$6.550316e-15$|$1.104750e-11$|$3.112093e-15$|$1024$|$0.25$|
|$128$|$8.326672e-15$|$8.065748e-11$|$2.220446e-16$|$4096$|$0.25$|
|$512$|$7.072120e-14$|$5.437144e-9$|$3.964173e-14$|$65536$|$0.25$|
|$1024$|$1.445787e-13$|$1.967658e-8$|$4.568075e-14$|$262144$|$0.25$|
|$4096$|$7.660539e-13$|$2.223743e-6$|$1.829228e-13$|$4194304$|$0.25$|
|$8192$|$9.765799e-13$|$2.502347e-5$|$5.187150e-14$|$16777216$|$0.25$|
Sign in with Wallet
Connect another wallet