| jupytext | kernelspec | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
FIZ228 - Numerical Analysis
Dr. Emre S. Tasci, Hacettepe University
+++
This lecture is heavily benefited from Steven Chapra's [Applied Numerical Methods with MATLAB: for Engineers & Scientists](https://www.mheducation.com/highered/product/applied-numerical-methods-matlab-engineers-scientists-chapra/M9780073397962.html).
+++
In our "Regression" lecture, we talked about fitting a function/model/curve that passes as close as possible to a given set of data points. But sometimes we find ourselves in a case where we have clearly defined points and would like to know about the behaviour between these points. This method of estimating the intermediate values while retaining the known values is known as interpolation.
The given data points act as the boundary conditions, usually fixed and thus very suitable for interpolation applcations. We are given a partial look on the system which corresponds to the initial conditions.
As was the case in regression, there is always a risk of overfitting alas in the interpolation, we must hit all the given data (within an error tolerance that is related to the measurements.
import numpy as np
import matplotlib.pyplot as plt
Revisiting our discussion in regression, we can always fit a set of
| 1 | 10 | 25 |
| 2 | 20 | 70 |
| 3 | 30 | 380 |
| 4 | 40 | 550 |
| 5 | 50 | 610 |
| 6 | 60 | 1220 |
| 7 | 70 | 830 |
| 8 | 80 | 1450 |
+++
We can use numpy's built in polynomial functions to fit (polyfit), evaluate (poly1d & polyval) and even construct polynomials from roots (poly) or vice versa (roots).
data = np.array([range(10,90,10),[25,70,380,550,610,1220,830,1450]]).T
x = data[:,0]
y = data[:,1]
p = np.polyfit(x,y,len(x)-1)
print(p)
xx = np.linspace(10,80,100)
yy = np.zeros(len(xx))
n = len(x)
for k in range(n):
yy += p[k]*xx**(n-k-1)
print("-"*70)
# we could as well had used poly1d function
# to functionalize the polynomial 8)
f = np.poly1d(p)
print(f)
print(f(x))
# another alternative is to use polyval:
g = np.polyval(p,x)
print(g)
plt.plot(xx,yy,"-b",x,y,"ok",xx,f(xx),"-r")
plt.show()
A 7th degree polynomial might be too much for our taste! 8) So let's begin with something simpler. For example, what would be the polynomial whose roots are -1 and 2?..
Since it has two roots, it must be a second order polynomial (quadratic):
q = np.poly([-1,2])
print(q)
To write it more neatly, we can construct a poly object using the coefficients:
qq = np.poly1d(q)
print(qq)
ladies and gents!
Let's go from reverse and find the roots of a given polynomial via roots:
x12 = np.roots(qq)
print(x12)
# feeding the coefficients would also work:
x12 = np.roots(q)
print(x12)
xx = np.linspace(-5,5,100)
yy = np.polyval(qq,xx)
plt.plot(xx,yy,"-r")
plt.xticks(np.arange(-5,6,1))
plt.yticks(np.arange(0,26,5))
plt.grid(True)
plt.show()
Now, let's take 3 points from our designated polynomial, and afterwards forget about the actual polynomial:
qc = [1, -1, 2]
qq = np.poly1d(qc)
x = [-3.1,0.7,4.3]
y = qq(x)
print(np.array([x,y]).T)
xx = np.linspace(-5,5,100)
yy = np.polyval(qq,xx)
plt.plot(xx,yy,color="gray")
plt.plot(x,y,"ko")
plt.xticks(np.arange(-5,6,1))
plt.yticks(np.arange(0,26,5))
plt.grid(True)
plt.show()
This is the part where we forget our polynomial information:
and thus, we are left with the data points to ponder on:
| x | -3.1 | 0.7 | 4.3 |
|---|---|---|---|
| y | 14.71 | 1.79 | 16.19 |
We have 3 data-points, so we'll assume a 2nd order polynomial in the form:
and we have:
x = np.array([-3.1, 0.7, 4.3])
print(x**2)
so, it is actually nothing but 3 equations with 3 unknowns! We know how to solve it! 8)
$$\begin{pmatrix} 9.61 & -3.1 & 1\ 0.49&0.7&1\ 18.49&4.3&1 \end{pmatrix} \begin{pmatrix}q_1\q_2\q_3\end{pmatrix}
\begin{pmatrix}14.71\1.79\16.19\end{pmatrix} $$
A = np.array([[9.61,-3.1,1],[0.49,0.7,1],[18.49,4.3,1]])
b = np.array([14.71,1.79,16.19])
q123 = np.linalg.solve(A,b)
print(q123)
... ta-taaaa!
+++
Newton polynomials incorporate known (given) points to go for the next one. For example, if two points are known, we pass a line through them and using similarity, we evaluate the value at the specified point.
Suppose we have
(don't forget that we don't know / can't see the real function but only know about the two data points $(x_1,x_2)$)
So, we'll reach
singling out
The indice "1" in $f_1$ denotes that this is a first-order (linear) interpolation.
+++
Estimate
-
$\ln1 = 0$ and$\ln6=1.791759$ -
$\ln1 = 0$ and$\ln4=1.386294$
(calculate the true percentage relative error (
# 1
ln1 = 0
ln6 = 1.791759
ln2_1 = ln1 + (ln6-ln1)/(6-1)*(2-1)
print(" ln2 ~ {:.6f}".format(ln2_1))
print("True value: {:.6f}".format(np.log(2)))
print(" E_t: {:.2f}%".format((np.log(2)-ln2_1)/np.log(2)*100))
# 2
ln1 = 0
ln4 = 1.386294
ln2_2 = ln1 + (ln4-ln1)/(4-1)*(2-1)
print(" ln2 ~ {:.6f}".format(ln2_2))
print("True value: {:.6f}".format(np.log(2)))
print(" E_t: {:.2f}%".format((np.log(2)-ln2_2)/np.log(2)*100))
Newton's polynomials are represented in increasing order as:
here, the
these coefficients are called as divided differences.
Here's a scheme of how to proceed forward in calculating these divided differences:
Let's put this into action in our next (evolved) example:
+++
Use the three of the previous given
$x_1 = 1,;f(x_1) = 0$ $x_2 = 4,;f(x_2)=1.386294$ $x_3 = 6,;f(x_3)=1.791759$
(calculate the true percentage relative error (
+++
x1 = 1
f_x1 = 0
x2 = 4
f_x2 = 1.386294
x3 = 6
f_x3 = 1.791759
b1 = f_x1
b2 = (f_x2 - f_x1)/(x2-x1)
c2 = (f_x3 - f_x2)/(x3-x2)
b3 = (c2 - b2)/(x3-x1)
x = 2
f_x = f_x1 + b2*(x-x1)+b3*(x-x1)*(x-x2)
print(" ln2 ~ {:.6f}".format(f_x))
print("True value: {:.6f}".format(np.log(2)))
print(" E_t: {:.2f}%".format((np.log(2)-f_x)/np.log(2)*100))
# Plotting the given data, our estimation for x=2
# and the actual function:
xx = np.linspace(1,7,300)
yy = np.log(xx)
yN = f_x1 + b2*(xx-x1)+b3*(xx-x1)*(xx-x2)
data_x = np.array([x1,x2,x3,x])
data_y = np.array([f_x1,f_x2,f_x3,f_x])
plt.plot(xx,yy,"-",color="lightgray")
plt.plot(xx,yN,"-",color="mistyrose")
plt.plot(data_x[0:3],data_y[0:3],"ob")
plt.plot(data_x[3],data_y[3],"or")
plt.show()
The advantage of this method is its flexibility in the face of the addition of more data points. So suppose that we made another measurement at the
x1 = 1
f_x1 = 0
x2 = 4
f_x2 = 1.386294
x3 = 6
f_x3 = 1.791759
b1 = f_x1
b2 = (f_x2 - f_x1)/(x2-x1)
c2 = (f_x3 - f_x2)/(x3-x2)
b3 = (c2 - b2)/(x3-x1)
# --- up to here, we had the same calculations,
# no need to recalculate them in actual run ---
x4 = 5
f_x4 = 1.609438
# we'll fill in the new coefficients:
d2 = (f_x4 - f_x3)/(x4-x3)
c3 = (d2 - c2)/(x4-x2)
b4 = (c3 - b3)/(x4-x1)
# and we just append the b4 coefficient in our approximation:
x = 2
f_x = f_x1 + b2*(x-x1)+b3*(x-x1)*(x-x2)\
+b4*(x-x1)*(x-x2)*(x-x3)
print(" ln2 ~ {:.6f}".format(f_x))
print("True value: {:.6f}".format(np.log(2)))
print(" E_t: {:.2f}%".format((np.log(2)-f_x)/np.log(2)*100))
# Plotting the given data, our estimation for x=2
# and the actual function:
xx = np.linspace(1,7,300)
yy = np.log(xx)
yN = f_x1 + b2*(xx-x1)+b3*(xx-x1)*(xx-x2)+b4*(xx-x1)*(xx-x2)*(xx-x3)
data_x = np.array([x1,x2,x3,x4,x])
data_y = np.array([f_x1,f_x2,f_x3,f_x4,f_x])
plt.plot(xx,yy,"-",color="black")
plt.plot(xx,yN,"-",color="red")
plt.plot(data_x[0:4],data_y[0:4],"ob")
plt.plot(data_x[4],data_y[4],"or")
plt.show()
Lagrange polynomials are yet another way to interpolate using the given data points where it uses the addition of the weighted data points.
(Figure: Chapra, 2012)
1st order Lagrange polynomial is given by:
where:
so:
2nd order Lagrange polynomial is given by:
where:
so:
and in general:
+++
Use a Lagrange interpolating polynomial of the first and second order to evaluate the density of unused motor oil at
x1 = 0
f_x1 = 3.85
x2 = 20
f_x2 = 0.8
x = 15
f1_x = (x-x2) / (x1-x2) * f_x1\
+(x-x1) / (x2-x1) * f_x2
print("1st order approximation: {:.6f}".format(f1_x))
x3 = 40
f_x3 = 0.212
f2_x = ((x-x2) * (x-x3)) / ((x1-x2) * (x1-x3)) * f_x1\
+((x-x1) * (x-x3)) / ((x2-x1) * (x2-x3)) * f_x2\
+((x-x1) * (x-x2)) / ((x3-x1) * (x3-x2)) * f_x3
print("2nd order approximation: {:.6f}".format(f2_x))
Use the given
$x_1 = 1,;f(x_1) = 0$ $x_2 = 4,;f(x_2)=1.386294$ $x_3 = 6,;f(x_3)=1.791759$
(calculate the true percentage relative error (
x1 = 1
f_x1 = 0
x2 = 4
f_x2 = 1.386294
x3 = 6
f_x3 = 1.791759
x = 2
f2_x = ((x-x2) * (x-x3)) / ((x1-x2) * (x1-x3)) * f_x1\
+((x-x1) * (x-x3)) / ((x2-x1) * (x2-x3)) * f_x2\
+((x-x1) * (x-x2)) / ((x3-x1) * (x3-x2)) * f_x3
print(" ln2 ~ {:.6f}".format(f2_x))
print("True value: {:.6f}".format(np.log(2)))
print(" E_t: {:.2f}%".format((np.log(2)-f2_x)/np.log(2)*100))
## Using the built in implementation:
from scipy.interpolate import lagrange
x = [1,4,6]
y = [0,1.386294, 1.791759]
p = lagrange(x,y)
print(p)
np.polyval(p,2)
Normally, we do interpolation to calculate
The first thing that comes to mind is to reinterpret
np.set_printoptions(precision=3)
x = np.arange(1,8)
y = 1/x
print(np.array([x,y]))
now suppose that we want to find the value of
print(np.array([y,x]))
p = np.polyfit(y,x,6)
pol = np.poly1d(p)
print(p)
print(pol)
np.polyval(p,0.3)
... or, we can use our heads, as usual! ;)
What we lacked in the above approximation is we failed to read the data. The data, with the given (x,y) order was evenly spaced. Therefore we could have gone for a 2nd order fit around
q = np.polyfit(x[1:4],y[1:4],2)
pol2 = np.poly1d(q)
print(pol2)
We are looking at a value of
and this is just a problem of finding the root!
x = np.roots([0.04167,-0.375,1.083-0.3])
print(x)
Continuing to using our heads, we remember that, we were looking for the solutions between 2 and 4, hence we take 3.293 -- considering that the true value being 3.333, we did a great job, by keeping to our human sides, not by being a number cruncher! 8)
+++
(not quite interpolation but still... ;)
When we have a differential equation and the values at boundaries are known, we can interpolate via the finite differences method to approximate the derivatives as the slopes of the lines connecting the data points.
Suppose that we want to calculate the derivative of a function at the
(Chapra, 2012)
a) Forward Finite Difference Approximation
Here, we have the slope as:
b) Backward Finite Difference Approximation
Considering the backward we can also write the slope as:
but the best of the both worlds would be to take both into account:
Centered Finite Difference Approximation
(Chapra, 2012)
+++
The figure shows a long, thin rod positioned between two walls that are held at constant temperatures. Heat flows through the rod as well as between the rod and the surrounding air. For the steady-state case, a differential equation based on heat conservation can be written for such a system as:
where
Given
+++
Solution
If we approximate the 2nd order derivative with the finite difference method, we will have for each of the marked points (
and the main differential equation takes the form of:
Expanding and re-arranging yields:
substituting the values of
Writing explicitly for each of the points, we have:
We already have
This is nothing but a linear equation set! So all of a sudden our 2nd order ODE has transformed into a set of 4 linear equations with 4 unknowns!
A1 = np.diag([2.04,2.04,2.04,2.04])
A2 = np.diag([-1,-1,-1],1)
A3 = np.diag([-1,-1,-1],-1)
A = A1+A2+A3
print(A)
b = np.array([40.8,0.8,0.8,200.8])
print(b)
T = np.linalg.solve(A,b)
print(T)
Btw, the analytical solution of the given ODE is:
Let's check how good did we do:
def fT(x):
return 73.4523*np.exp(0.1*x)-53.4523*np.exp(-0.1*x)+20
xx = np.linspace(0,12,100)
TT = fT(xx)
plt.plot(xx,TT,"-",color="black")
plt.plot(np.arange(2,9,2),T,"ro")
plt.plot([0,10],[40,200],"ko")
plt.show()
- Heavily benefitted from Chapra's "Applied Numerical Methods with MATLAB for Engineers and Scientists" (2012) along with dear Dr. Eda Çelik Akdur's "KMU 231 - Numerical Analysis" lecture notes.
- The divided differences scheme's
$\LaTeX$ code was adapted from the corresponding Wikipedia entry.