---
title: 11-3. Chebfun. Non-standard boundary conditions
tags: Final reports
---

**Reference**
1. Chebfun - Nonstandard 'boundary' conditions, [NonstandardBCs](https://www.chebfun.org/examples/ode-linear/NonstandardBCs.html)

---

### Problem 1
Solve the ODE
$$u''+x^2u=1$$
on the domain $[-1,1]$ subject to the conditions
$$u(-1)=1; \text{ The average of u over } [-1,1] \text{ is } 0$$

#### solution
We solve the system
$$Mu=r$$
where
$$
M=\begin{pmatrix}
1 & 0 & \dots & 0 & 0\\
& & CCW\\
& & P\times(D^2+diag(x_i^2))\\
\end{pmatrix}_{N \times N}
\quad
u=\begin{pmatrix}
u(x_0)\\
u(x_1)\\
u(x_2)\\
\vdots\\
u(x_{N-1})\\
\end{pmatrix}_{N \times 1}
\quad
r=\begin{pmatrix}
1\\
0\\
1\\
\vdots\\
1\\
\end{pmatrix}_{N \times 1}
$$
$CCW_{1 \times N}$ is the Clenshaw-Curtis quadrature for $N$ Chebyshev points of the second kind on $[-1,1]$.
$D_{N \times N}$ is the differentiation matrix of $N$ Chebyshev points of the second kind on $[-1,1]$.
$P_{(N-2) \times N}$ is the Barycentric resampling matrix from $N$ Chebyshev points of the second kind on $[-1,1]$ to $N-2$ Chebyshev points of the first kind on $[-1,1]$.

In particular,
$N=19$
$X: N$ Chebyshev points of the second kind on $[-1,1]$
$Y: N-2$ Chebyshev points of the first kind on $[-1,1]$
$$
X=\begin{pmatrix}
-cos(\frac{0}{N-1})\\
-cos(\frac{\pi}{N-1})\\
-cos(\frac{2\pi}{N-1})\\
\vdots\\
-cos(\frac{(N-1)\pi}{N-1})\\
\end{pmatrix}_{N \times 1}
\quad
Y=\begin{pmatrix}
-cos(\frac{(2 \times 0+1)\pi}{2 \times (N-3)+2})\\
-cos(\frac{(2 \times 1+1)\pi}{2 \times (N-3)+2})\\
-cos(\frac{(2 \times 2+1)\pi}{2 \times (N-3)+2})\\
\vdots\\
-cos(\frac{(2 \times (N-3)+1)\pi}{2 \times (N-3)+2})\\
\end{pmatrix}_{(N-2) \times 1}
$$
$$
[P]_{j,k} = \left\{\begin{array}{ll}
\frac{\lambda_k}{y_j-x_k}(\sum_{l=0}^{N-1}{\frac{\lambda_l}{y_j-x_l}})^{-1} & \mbox{if $y_j \neq x_k$} \\ \mbox
1 & \mbox{if $y_j = x_k$}
\end{array} \right.
\quad
\lambda_k = \left\{\begin{array}{ll}
\frac{1}{2} & \mbox{if $k = N-1$} \\ \mbox
(-1)^{N-1-k} & \mbox{if $k=1, \dots,N-2$} \\
\frac{(-1)^{N-1}}{2} & \mbox{if $k=0$}
\end{array} \right.
$$

#### Result:

Chebyshev coefficients:

0.0829307517870041
0.682831733730336
0.249115851228505
-0.0108499779455154
-0.0013955717700003
-0.00213364601348834
-0.000520017396039883
1.69844980867949e-005
1.68997381282981e-006
1.85206724078673e-006
3.61366047563883e-007
-9.83598570380614e-009
-8.27232621137942e-010
-7.4200335208669e-010
-1.2413734238079e-010
2.95719159036799e-012
2.1896065207885e-013
1.68951272725179e-013
2.51218798849903e-014

Reference:
[CCW code](https://github.com/chebfun/chebfun/blob/master/%40chebtech2/quadwts.m)

### Problem 2
Solve the ODE
$$u''+x^2u=1$$
on the domain $[-1,1]$ subject to the conditions
$$u(-1)=1; \quad u(0)=0.5$$

#### solution
We solve the system
$$Mu=r$$
$$
M=\begin{pmatrix}
1 & 0 & \dots & 0 & 0 & 0 & \dots & 0 & 0\\
& & & & P\times(D^2+diag(x_i^2))\\
0 & 0 & \dots & 0 & 1 & 0 & \dots & 0 & 0\\
\end{pmatrix}_{N \times N}
\quad
u=\begin{pmatrix}
u(x_0)\\
u(x_1)\\
u(x_2)\\
\vdots\\
u(x_{N-1})\\
\end{pmatrix}_{N \times 1}
\quad
r=\begin{pmatrix}
1\\
1\\
\vdots\\
1\\
0.5\\
\end{pmatrix}_{N \times 1}
$$
$N$ is odd.
$D_{N \times N}$ is the differentiation matrix of $N$ Chebyshev points of the second kind on $[-1,1]$.
$P_{(N-2) \times N}$ is the Barycentric resampling matrix from $N$ Chebyshev points of the second kind on $[-1,1]$ to $N-2$ Chebyshev points of the first kind on $[-1,1]$.

#### Results:

Chebyshev coefficients:

0.729414210616675
0.0585218718079899
0.221753098829177
-0.000929893833402149
-0.00812877466199519
-0.000182863438116363
-0.000458780562370179
1.45565088826008e-006
9.20442339104293e-006
1.58730820882309e-007
3.17836254702674e-007
-8.4298989298481e-010
-4.38514892299303e-009
-6.35930617699209e-011
-1.09015782512323e-010
2.53485118697741e-013
1.14561446748793e-012
1.43773881688958e-014
2.20070875103475e-014

code 1 (C++)

//compile with g++ NonStandardBC.cpp -L.\. -lfftw3 -o NBC
#include <iostream>
#include <cmath>
#include <fstream>
#include <complex>
#include <vector>

#include "fftw3.h"

#define REAL 0
#define IMAG 1
#define PI 3.141592653589793

using namespace std;
void norm_real( fftw_complex c[], int N);  //normalize and set imaginary part to 0

int main(){

    ofstream outfile("./output.txt");
    outfile.precision(15);

    // Set number of points
    int N = 19;

    // Compute Differentiation Matrix D
    vector<double> x(N,0);    // N chebyshev points of the second kind on [-1,1]
    for (int i=0;i<N;i++) { x[i] = -cos((PI*i) / (N-1)); }
	
    vector<double> c(N,1);
    c[0] = 2; c[N-1] = 2;
    for (int i=0;i<N;i++) { c[i] = c[i]*pow((-1.0),i); }
	
    double X[N][N],Y[N][N];
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		X[j][i] = x[j];
    		Y[i][j] = x[j];    //Y = X transpose
    	}
    }
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		X[i][j] = X[i][j] - Y[i][j];   //dX = X-X'
    	}
    }
	
    double C[N][N],D[N][N];	
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		C[i][j] = c[i]*(1.0/c[j]);
    	}
    }
    for(int i=0;i<N;i++) { X[i][i] = X[i][i] + 1.0; }
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		D[i][j] = C[i][j]/X[i][j];
    	}
    }
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		Y[i][j] = D[j][i];   //Y = D transpose
    	}
    }

    vector<double> sum(N,0);
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    	    sum[i] = sum[i] + Y[j][i];   //Y = D transpose
	    }
    }
    for(int i=0;i<N;i++) { D[i][i] = D[i][i] - sum[i]; }
    
    // Compute lambda: coeff of Barycentric interpolation
    vector<double> lambda(N,1);
    lambda[0] = 0.5; lambda[N-1] = 0.5;
    for (int i=0;i<N;i++){ lambda[N-1-i] = lambda[N-1-i]*pow((-1.0),i); }

    // prepare for ifft
    fftw_complex *in, *out;
    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N-1));
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N-1));

    // Compute CCW
    for(int i=1; i<=floor((N-1)/2); i++)
    {
       in[i][REAL] = 2/(1-pow(2*i,2));
       in[i][IMAG] = 0;
    }
    for(int i=floor((N)/2)-1; i>=1; i--)
    {
        in[N-1-i][REAL] = in[i][REAL];
        in[N-1-i][IMAG] = 0;
    }
    in[0][REAL] = 2;
    in[0][IMAG] = 0;
    

    fftw_plan plan1 = fftw_plan_dft_1d(N-1, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(plan1);
    norm_real(out, N-1);

    vector<double> CCW(N,0);
    for(int i=0; i<N-1; i++)
    {
        CCW[i] = out[i][REAL];
    }
    CCW[N-1] = CCW[0]/2;
    CCW[0] = CCW[0]/2;
    

    // N-2 chebyshev points of the first kind on [-1,1]
    vector<double> xx(N,0);    
    for (int i=0;i<N-2;i++){
    	xx[i] = -cos(PI*(2*i+1) / (2*(N-3)+2));
    }

    // Compute P: Barycentric resampling matrix
    double P[N-2][N];
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            P[i][j] = xx[i]-x[j];
        }
    }
    vector<double> rowsum(N-2,0);
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            P[i][j] = lambda[j]/P[i][j];
            rowsum[i] += P[i][j];
        }
        rowsum[i] = 1/rowsum[i];
    }
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            P[i][j] = P[i][j]*rowsum[i];
            if(isnan(P[i][j])) P[i][j] = 1;
        }
    }

    // Compute PA = P*(D^2 + diag[x^2]) and M
    double A[N][N]={0};
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            A[i][j] = 0;
        }
        A[i][i] = x[i]*x[i];
    }
    
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            for(int k=0; k<N; k++){
                A[i][j] += D[i][k]*D[k][j];
            }
        }
    }

    double PA[N-2][N]={0};
    double M[N][N]={0};
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            for(int k=0; k<N; k++){
                PA[i][j] += P[i][k]*A[k][j];
            }
            M[i+2][j] = PA[i][j];
        }
    }
    M[0][0] = 1;
    for(int i=1; i<2; i++){
        for(int j=0; j<N; j++){
            M[i][j] = CCW[j];
        }
    }

    // Set rhs
    vector<double> rhs(N,1);
    rhs[1]=0;

    // solve u
    // Gaussian elimination
    for(int i=0; i<N-1; i++){
        for(int j=i+1; j<N; j++){
            for(int k=i+1; k<N; k++){
                M[j][k] = M[j][k] - (M[j][i]/M[i][i])*M[i][k];
            }
            rhs[j] = rhs[j] - (M[j][i]/M[i][i])*rhs[i];
        }
    }

    // backward substitution
    vector<double> u(N,0);
    double s;
    u[N-1] = rhs[N-1]/M[N-1][N-1];
    for(int i=N-2; i>=0; i--){
        s = 0;
        for(int j=i+1; j<N; j++){
            s += M[i][j]*u[j];
        }
        u[i] = (rhs[i]-s)/M[i][i];
    }

    // output u
    outfile << "[" ;
    for(int i=0; i<N-1; i++)
    {
        outfile << u[i] << ", ";
    }
    outfile << u[N-1] << "]";
    
    fftw_destroy_plan(plan1);
    fftw_free(in);
    fftw_free(out);
    
    // Compute chebyshev coefficients
    fftw_complex *in2, *out2;
    in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2*N-2));
    out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2*N-2));
    fftw_plan plan2 = fftw_plan_dft_1d(2*N-2, in2, out2, FFTW_FORWARD, FFTW_ESTIMATE);
    
    for(int i=0; i<N; i++){
        in2[i][REAL] = u[i];
        in2[i][IMAG] = 0;
    }
    for(int i=N; i<(2*N-2); i++){
        in2[i][REAL] = u[2*N-2-i];
        in2[i][IMAG] = 0;
    }

    fftw_execute(plan2);
    norm_real(out2, 2*N-2);
    vector<double> coeffs(N,0);
    coeffs[0] = out2[0][REAL];
    for(int i=1; i<N-1; i++){
        coeffs[i] = 2*out2[i][REAL];
    }
    coeffs[N-1] = out2[N-1][REAL];

    // output coefficients
    outfile << endl << "chebyshev coefficients:" << endl;
    for(int i=0; i<N; i++){
        outfile << coeffs[i] << endl;
    }

    fftw_destroy_plan(plan2);
    fftw_free(in2);
    fftw_free(out2);

    return 0;
}

void norm_real( fftw_complex c[], int N)
{
    for(int i=0; i<N; i++)
    {
        c[i][REAL] = c[i][REAL]/N;
        c[i][IMAG] = 0;
    }
}

code 2(C++)

//compile with g++ NonStandardBC_EX2.cpp -L.\. -lfftw3 -o NBC2
#include <iostream>
#include <cmath>
#include <fstream>
#include <complex>
#include <vector>

#include "fftw3.h"

#define REAL 0
#define IMAG 1
#define PI 3.141592653589793

using namespace std;
void norm_real( fftw_complex c[], int N);  //normalize and set imaginary part to 0

int main(){

    ofstream outfile("./output.txt");
    outfile.precision(15);

    // Set number of points
    int N = 19;

    // Compute Differentiation Matrix D
    vector<double> x(N,0);    // N chebyshev points of the second kind on [-1,1]
    for (int i=0;i<N;i++) { x[i] = -cos((PI*i) / (N-1)); }
	
    vector<double> c(N,1);
    c[0] = 2; c[N-1] = 2;
    for (int i=0;i<N;i++) { c[i] = c[i]*pow((-1.0),i); }
	
    double X[N][N],Y[N][N];
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		X[j][i] = x[j];
    		Y[i][j] = x[j];    //Y = X transpose
    	}
    }
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		X[i][j] = X[i][j] - Y[i][j];   //dX = X-X'
    	}
    }
	
    double C[N][N],D[N][N];	
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		C[i][j] = c[i]*(1.0/c[j]);
    	}
    }
    for(int i=0;i<N;i++) { X[i][i] = X[i][i] + 1.0; }
	for(int i=0;i<N;i++){
            for(int j=0;j<N;j++){
		D[i][j] = C[i][j]/X[i][j];
	    }
    }
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		Y[i][j] = D[j][i];   //Y = D transpose
    	}
    }

    vector<double> sum(N,0);
    for(int i=0;i<N;i++){
    	for(int j=0;j<N;j++){
    		sum[i] = sum[i] + Y[j][i];   //Y = D transpose
    	}
    }
    for(int i=0;i<N;i++) { D[i][i] = D[i][i] - sum[i]; }
    
    // Compute lambda: coeff of Barycentric interpolation
    vector<double> lambda(N,1);
    lambda[0] = 0.5; lambda[N-1] = 0.5;
    for (int i=0;i<N;i++){ lambda[N-1-i] = lambda[N-1-i]*pow((-1.0),i); }

    // N-2 chebyshev points of the first kind on [-1,1]
    vector<double> xx(N,0);    
    for (int i=0;i<N-2;i++){
    	xx[i] = -cos(PI*(2*i+1) / (2*(N-3)+2));
    }

    // Compute P: Barycentric resampling matrix
    double P[N-2][N];
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            P[i][j] = xx[i]-x[j];
        }
    }
    vector<double> rowsum(N-2,0);
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            P[i][j] = lambda[j]/P[i][j];
            rowsum[i] += P[i][j];
        }
        rowsum[i] = 1/rowsum[i];
    }
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            P[i][j] = P[i][j]*rowsum[i];
            if(isnan(P[i][j])) P[i][j] = 1;
        }
    }

    // Compute PA = P*(D^2 + diag[x^2]) and M
    double A[N][N]={0};
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            A[i][j] = 0;
        }
        A[i][i] = x[i]*x[i];
    }
    
    for(int i=0; i<N; i++){
        for(int j=0; j<N; j++){
            for(int k=0; k<N; k++){
                A[i][j] += D[i][k]*D[k][j];
            }
        }
    }

    double PA[N-2][N]={0};
    double M[N][N]={0};
    for(int i=0; i<N-2; i++){
        for(int j=0; j<N; j++){
            for(int k=0; k<N; k++){
                PA[i][j] += P[i][k]*A[k][j];
            }
            M[i+1][j] = PA[i][j];
        }
    }
    M[0][0] = 1;
    M[N-1][N/2] = 1;

    // Set rhs
    vector<double> rhs(N,1);
    rhs[N-1]=0.5;

    // solve u
    // Gaussian elimination
    for(int i=0; i<N-1; i++){
        for(int j=i+1; j<N; j++){
            for(int k=i+1; k<N; k++){
                M[j][k] = M[j][k] - (M[j][i]/M[i][i])*M[i][k];
            }
            rhs[j] = rhs[j] - (M[j][i]/M[i][i])*rhs[i];
        }
    }

    // backward substitution
    vector<double> u(N,0);
    double s;
    u[N-1] = rhs[N-1]/M[N-1][N-1];
    for(int i=N-2; i>=0; i--){
        s = 0;
        for(int j=i+1; j<N; j++){
            s += M[i][j]*u[j];
        }
        u[i] = (rhs[i]-s)/M[i][i];
    }

    // output u
    outfile << "[" ;
    for(int i=0; i<N-1; i++)
    {
        outfile << u[i] << ", ";
    }
    outfile << u[N-1] << "]" << endl;

    // Compute chebyshev coefficients
    fftw_complex *in2, *out2;
    in2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2*N-2));
    out2 = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (2*N-2));
    fftw_plan plan2 = fftw_plan_dft_1d(2*N-2, in2, out2, FFTW_FORWARD, FFTW_ESTIMATE);
    
    for(int i=0; i<N; i++){
        in2[i][REAL] = u[i];
        in2[i][IMAG] = 0;
    }
    for(int i=N; i<(2*N-2); i++){
        in2[i][REAL] = u[2*N-2-i];
        in2[i][IMAG] = 0;
    }

    fftw_execute(plan2);
    norm_real(out2, 2*N-2);
    vector<double> coeffs(N,0);
    coeffs[0] = out2[0][REAL];
    for(int i=1; i<N-1; i++){
        coeffs[i] = 2*out2[i][REAL];
    }
    coeffs[N-1] = out2[N-1][REAL];

    // output coefficients
    outfile << endl << "chebyshev coefficients:" << endl;
    for(int i=0; i<N; i++){
        outfile << coeffs[i] << endl;
    }

    fftw_destroy_plan(plan2);
    fftw_free(in2);
    fftw_free(out2);

    return 0;
}

void norm_real( fftw_complex c[], int N)
{
    for(int i=0; i<N; i++)
    {
        c[i][REAL] = c[i][REAL]/N;
        c[i][IMAG] = 0;
    }
}