---
# System prepended metadata

title: 11-6. Spectral integration and two point BVP
tags: [Final reports]

---

---
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$|



