---
title: 9-3. Chebyshev approximations - interpolation
tags: Spectral method
GA: G-RYTKXJR91W
---

**Reference**
1. [Approximation theory and approximation practice](https://people.maths.ox.ac.uk/trefethen/ATAP/ATAPfirst6chapters.pdf) by *Lloyd N. Trefethen*

---

## Exercise 1

Consider
$$
f(x) = e^x \sin(300 x), \quad -1\le x\le 1.
$$

1. Use "enough" grid points to approximation the function and its derivative.

2. Write a program to interpolate f'(x) at given array of points. That is, the input of your program is a vector $\vec{x}$ and the output is the vector $f'(\vec{x})$.

* Randomly choose $10^3$, $10^4$ and $10^5$ points on $[-1, 1]$, how long does it takes for the interpolation?

---

###CODE(PYTHON):

import math
import numpy as np
import random as random
import matplotlib.pyplot as plt

#grid number
n=330

X=np.cos(np.linspace(0,n,n)*np.pi/n)
F=np.exp(X)*np.sin(300*X)

#Barycentric interpolation in Chebyshev points
def Baryc(x,n,X,F):
    a=((-1)**0*F[0]/(x-X[0])+(-1)**n*F[0]/(x-X[n-1]))*0.5
    for j in range(1,n-1):
        a=a+(-1)**j*F[j]/(x-X[j])
    b=((-1)**0/(x-X[0])+(-1)**n/(x-X[n-1]))*0.5
    for j in range(1,n-1):
        b=b+(-1)**j/(x-X[j])
    return a/b

#First order differential  Barycentric interpolation in Chebyshev points   
def D_Baryc(x,n,X,F):
    a=(-(-1)**0*F[0]/(x-X[0])**2-(-1)**n*F[0]/(x-X[n-1])**2)*0.5
    for j in range(1,n-1):
        a=a-(-1)**j*F[j]/(x-X[j])**2
    b=((-1)**0/(x-X[0])+(-1)**n/(x-X[n-1]))*0.5
    for j in range(1,n-1):
        b=b+(-1)**j/(x-X[j])
    A=a/b

#solve f(x)
Test_10000=2*np.random.rand(1,10000)-1
Test_10000=np.sort(Test_10000[0])
FTest_10000=Baryc(Test_10000,n,X,F)
plt.plot(Test_10000,FTest_10000)

#solve f'(x)
Test_10000=2*np.random.rand(1,10000)-1
Test_10000=np.sort(Test_10000[0])
FTest_10000=D_Baryc(Test_10000,n,X,F)
plt.plot(Test_10000,FTest_10000)

###############################################################
#Chebyshev coefficients

# generate grid points
N = 380;
X = np.cos(np.arange(0,N)*np.pi/N)

# evaluate at grid points
f = np.exp(X)*np.sin(300*X)

# extend the function and take scaled FFT

res=f[1:N-1]
res=res[::-1]
f_extend=np.hstack((f,res))

f_hat = fft(f_extend)/(2*N)

# obtain the Chebyshev coefficients
a = np.zeros((1,N+1))
a[0][0] = np.real(f_hat[0])
a[0][1:N-1] = 2*np.real(f_hat[1:N-1])
a[0][N] = np.real(f_hat[N])
plt.plot(np.log10(abs(a[0])))

### ANS of Q1:

For this function,more than 340 grid points is enough.

f(x):
f'(x):

Error estimation:
$log (a)$:
$log_{10} (a)$:

### ANS of Q2:

#### Randomly choose 100 points
f(x):
f'(x):cost 0.00800sec
#### Randomly choose 1000 points
f(x):
f'(x):cost 0.01400sec
#### Randomly choose 10000 points
f(x):
f'(x): cost 0.05803sec
f'(x) on the range[-0.99,1]

=========================================================
(閔中)

(Matlab code)
function Barycentric

N = 1000;
x = cos((2*(N:-1:1)-1)*pi/(2*N))';
x_plot = linspace(min(x-sqrt(5)/10000),max(x+sqrt(5)/10000),2*N)';

f = @(x) exp(x).*sin(300*x);
df = @(x) exp(x).*(sin(300*x)+300*cos(300*x));

p = chebf(x,x_plot,f);
dp = chebdf(x,x_plot,df);
errorf = abs(p-f(x_plot));
errordf = abs(dp-df(x_plot));

figure(1)
subplot(3,1,1)
plot(x_plot,p,'r') ; hold on;
axis([-1 1 -3 3])
title('f(x)')
ylabel('p(x)')
subplot(3,1,2)
plot(x_plot,f(x_plot),'k'); hold on;
axis([-1 1 -3 3])
ylabel('f(x)')
subplot(3,1,3)
semilogy(x_plot,errorf,'b'); hold off;
axis([-1 1 10^-5 1])
ylabel('|f(x)-p(x)|')
xlabel('x') 

figure(2)
subplot(3,1,1)
plot(x_plot,dp,'r') ; hold on;
axis([-1 1 -1000 1000])
title('df/dx')
ylabel('dp/dx')
subplot(3,1,2)
plot(x_plot,df(x_plot),'k'); hold on;
axis([-1 1 -1000 1000])
ylabel('df/dx')
subplot(3,1,3)
semilogy(x_plot,errordf,'b'); hold off;
axis([-1 1 10^-5 1])
ylabel('|df/dx-dp/dx|')
xlabel('x') 

tic
x2 = -1+2*rand(100,1);
x2 = sort(x2);
p2 = chebf(x,x2,f);
dp2 = chebdf(x,x2,df);
toc
tic
x3 = -1+2*rand(1000,1);
x3 = sort(x3);
p3 = chebf(x,x3,f);
dp3 = chebdf(x,x3,df);
toc
tic
x4 = -1+2*rand(10000,1);
x4 = sort(x4);
p4 = chebf(x,x4,f);
dp4 = chebdf(x,x4,df);
toc
tic
x5 = -1+2*rand(100000,1);
x5 = sort(x5);
p5 = chebf(x,x5,f);
dp5 = chebdf(x,x5,df);
toc

figure(3)
subplot(4,1,1)
plot(x2,p2,'r') ; hold on;
axis([-1 1 -3 3])
title('p(x)')
ylabel('p(x) at n = 10^2')
subplot(4,1,2)
plot(x3,p3,'b') ; hold on;
axis([-1 1 -3 3])
ylabel('p(x) at n = 10^3')
subplot(4,1,3)
plot(x4,p4,'g'); hold on;
axis([-1 1 -3 3])
ylabel('p(x) at n = 10^4')
subplot(4,1,4)
plot(x5,p5,'k'); hold off;
axis([-1 1 -3 3])
ylabel('p(x) at n = 10^5')
xlabel('x') 

figure(4)
subplot(4,1,1)
plot(x2,dp2,'r') ; hold on;
axis([-1 1 -1000 1000])
title('dp/dx')
ylabel('dp/dx at n = 10^2')
subplot(4,1,2)
plot(x3,dp3,'b') ; hold on;
axis([-1 1 -1000 1000])
ylabel('dp/dx at n = 10^3')
subplot(4,1,3)
plot(x4,dp4,'g'); hold on;
axis([-1 1 -1000 1000])
ylabel('dp/dx at n = 10^4')
subplot(4,1,4)
plot(x5,dp5,'k'); hold off;
axis([-1 1 -1000 1000])
ylabel('dp/dx at n = 10^5')
xlabel('x') 

function p = chebf(x,x_plot,f)
N = length(x);
p = zeros(length(x_plot),1);
for m = 1:length(x_plot)
p(m)  = sum((-1).^(1:N)'.*f(x)./(x_plot(m)-x))/sum((-1).^(1:N)'./(x_plot(m)-x));
end

function dp = chebdf(x,x_plot,df)
N = length(x);
dp = zeros(length(x_plot),1);
for m = 1:length(x_plot)
dp(m)  = sum((-1).^(1:N)'.*df(x)./(x_plot(m)-x))/sum((-1).^(1:N)'./(x_plot(m)-x));
end

## Exercise 2

Consider
$$
g(x) = \sin(e^{10 x}), \quad -1\le x\le 1.
$$

1. Use "enough" grid points to approximation the function, denoted by $p(x)$. (It should be more than $10^4$)

2. Write a program to interpolate p(x) at $1000$ random points. That is, the input of your program is a vector $\vec{x}$ and the output is the vector $p(\vec{x})$.

* How long does it takes to evluate $p(\vec{x})$?
* How long does it takes to evluate $g(\vec{x})$?
* Draw a plot for $|p(\vec{x}) - g(\vec{x})|$

---
### Answer:
A1. p(x) cost 113.978306 sec

A2. g(x) cost 0.005985 sec

A3.


error estimation

Python code

import numpy as np
import matplotlib.pyplot as plt

#在-1到1之間取 N 個 Chebyshev points 和對應的值
N=42800
x=np.cos((np.linspace(0,N,N+1)*np.pi)/N)
f=np.sin(np.exp(10*x))

#Barycentric_interpolation
def Barycentric_interpolation(x0,x,f,N):
    frac=0
    deno=0
    for j in range(1,N):
        frac=frac+((-1)**j*f[j])/(x0-x[j])
        deno=deno+((-1)**j)/(x0-x[j])
    
    frac=frac+0.5*(((-1)**0*f[0])/(x0-x[0]))+0.5*(((-1)**(N)*f[N])/(x0-x[N]))
    deno=deno+0.5*(((-1)**0)/(x0-x[0]))+0.5*(((-1)**(N))/(x0-x[N]))
    return frac/deno                                        

a=f.tolist()[1:N+1]
a.reverse()
A=f.tolist()+a
f_extend=np.array(A)

f_hat=np.real(np.fft.fft(f_extend))

a =np.zeros(N+1);
a[0]= f_hat[0];
a[1:N] = 2*f_hat[1:N];
a[N] = f_hat[N];

plt.plot(np.log10(abs(a)))
plt.savefig('the number of coefficient.png')

import random
import time

P=[]
for i in range(1000):
    P.append(random.uniform(-1, 1))
P=np.array(P)
P.sort()

tStart = time.time()

Res_num=[]
for i in range(len(P)):
    Res_num.append(Barycentric_interpolation(P[i],x,f,N))
Res_num=np.array(Res_num)

tEnd = time.time()
print ("p(x) cost %f sec" % (tEnd - tStart))
print (tEnd - tStart)

tStart = time.time()

Res_true=[]
for i in range(len(P)):
    Res_true.append(np.sin(np.exp(10*P[i])))
Res_true=np.array(Res_true)

tEnd = time.time()
print ("g(x) cost %f sec" % (tEnd - tStart))
print (tEnd - tStart)

import matplotlib.pyplot as plt
plt.plot(P,abs(Res_true-Res_num))
plt.ylabel('|p(x)-g(x)|')
plt.savefig('P-g.png')