# Lecture 8: Interpolation¶

Date: 09/26/2017, Tuessday

```
In [1]:
```

```
format compact
```

## A simple example¶

Make 3 data points

```
In [2]:
```

```
xp = [-pi/2, 0, pi/2];
yp = [-1, 0, 1];
```

The two functions below both go through the 3 data points.

```
In [3]:
```

```
f = @(x) sin(x);
g = @(x) 2*x/pi;
```

There are ways to plot a function symbolically/analytically (for example), but those methods have a lot of limitations and you don’t have detailed controls on them.

So we stick to the most standard way of plotting: evaluate the function value on a lot of points to make the line look smooth.

```
In [4]:
```

```
x = linspace(-2,2,1e3); % a 1000 grid points from -2 to 2
```

Always remember to surpress the output (by `;`

) when defining this
kind of large array. MATLAB WILL print all the elements in an
array/matrix no matter how large it is. Sometimes your program will die
just because it wants to print \(10^{10}\) numbers.

Now we can plot the functions and data points.

```
In [5]:
```

```
%plot --size 300,200
hold on
plot(xp,yp,'o')
plot(x,f(x))
plot(x,g(x))
legend('data','sin(x)','2x/\pi','Location','northwest')
```

## Polynomial interpolation¶

```
In [6]:
```

```
% make some data points
n = 5;
px = [0 1 2 3 4];
py = [1 2 2 6 9];
```

n data point can be precisely fitted by a n-1 degree polynomial \(p=c_1+c_2x+...+c_nx^{n-1}\).

The coffeicients \(c=[c_1,c_2,...c_n]^T\) satisfy the equation

where \(y=[y_1,y_2,...,y_n]^T\) is the data points you want to fit, and V is the vandermode matrix only containing powers of \(x_k\), the data points.

```
In [7]:
```

```
% calculate the vander matrix by loop
V = zeros(n);
for j=1:5
V(:,j) = px.^(j-1);
end
V
```

```
V =
1 0 0 0 0
1 1 1 1 1
1 2 4 8 16
1 3 9 27 81
1 4 16 64 256
```

The built-in `vander`

is flipped left-right. Both forms are correct as
long as your algorithm is consistent with the matrix.

```
In [8]:
```

```
vander(px)
```

```
ans =
0 0 0 0 1
1 1 1 1 1
16 8 4 2 1
81 27 9 3 1
256 64 16 4 1
```

Now we can solve \(Vc=y\) by \(c=V^{-1}y\). In MATLAB backslash
form it is `c=V\y`

. Note that the actual code never computes
\(V^{-1}\), but use something like LU factorization/Gaussian
elimination to solve the system. (it is almost always a bad idea to
compute the inverse of a matrix). But thinking about \(V^{-1}y\)
helps you to remember the order of V and y in the command (e.g. is it
`V\y`

or `y\V`

?).

The code below throws an error because y is a row vector.

```
In [9]:
```

```
c = V\py
```

```
Error using \
Matrix dimensions must agree.
```

You need a column vector on the right side, just like how you write the equation mathematically.

```
In [10]:
```

```
c = V\py'
```

```
c =
1.0000
5.6667
-7.5833
3.3333
-0.4167
```

Now we have the coefficients \(c\), we can write the polynomial \(p=c_1+c_2x+...+c_nx^{n-1}\) as a MATLAB function.

Here’s a naive way to evaluate the polynomial. You should write a loop instead. Also consider Horner’s method to achieve optimal performance.

```
In [11]:
```

```
my_ploy = @(c,x) c(1) + c(2)*x + c(3)*x.^2 + c(4)*x.^3 + c(5)*x.^4 ;
```

```
In [12]:
```

```
% evaluate the function on a lot of data points
xf = linspace(-1,5,100);
yf = my_ploy(c,xf);
```

```
In [13]:
```

```
%plot --size 400,200
hold on
plot(px,py,'o')
plot(xf,yf)
legend('data','interpolation','Location','southeast')
```