Problem Set on Interpolation Using IPython
Q12
Assume $f(x)=x^2+3x+2,x\in[0,1]$, please calculate the greatest square approximation polynomial of $f(x)$ when $\rho (x)=1,\Phi=span\{1,x\}$. What if we take $\Phi=span\{1,x,x^2\}$?
Solution:
$$d_0=\int^1_0f(x)dx$$
%pylab inline
%config InlineBackend.figure_format = 'retina'
from sympy import *
init_printing()
x = Symbol("x")
integrate(x**2+3*x+2,(x,0,1))
Populating the interactive namespace from numpy and matplotlib
$$\frac{23}{6}$$
integrate(x*(x**2+3*x+2),(x,0,1))
$$\frac{9}{4}$$
A=Matrix([ [1, Rational(1,2) ],[Rational(1,2),Rational(1,3)] ])
A
$$\left[\begin{matrix}1 & \frac{1}{2}\\\frac{1}{2} & \frac{1}{3}\end{matrix}\right]$$
A.LUsolve(Matrix(2,1,[Rational(23,6),Rational(9,4)]))
$$\left[\begin{matrix}\frac{11}{6}\\4\end{matrix}\right]$$
thus the polynomial is $S_1^*(x)=\frac{11}{6}+4x$
for $\text{span}\{1,x,x^2\}$, similarily:
d3=integrate(x**2*(x**2+3*x+2),(x,0,1))
A=Matrix([ [1, Rational(1,2), Rational(1,3) ],[Rational(1,2),Rational(1,3),Rational(1,4)],[Rational(1,3),Rational(1,4),Rational(1,5)] ])
A.LUsolve(Matrix(3,1,[Rational(23,6),Rational(9,4),d3]))
$$\left[\begin{matrix}2\\3\\1\end{matrix}\right]$$
We can see that the result is exactly the same as the original $f(x)$.
Q14
Calculate the greatest square approx. polynomial of $f(x)$ on $\Phi=\text{span}\{1,x\}$ for:
- $f(x)=\frac{1}{x},[1,3]$
- $f(x)=e^x,[0,1]$
- $f(x)=\text{cos}\pi x,[0,1]$
- $f(x)=\text{ln}x,[1,2]$
def linearapprox(f,x0,x1):
d0=integrate(f,(x,x0,x1))
d1=integrate(x*(f),(x,x0,x1))
A=Matrix([ [1, Rational(1,2) ],[Rational(1,2),Rational(1,3)] ])
return A.LUsolve(Matrix(2,1,[d0,d1]))
linearapprox(1/x,1,3)
$$\left[\begin{matrix}-12 + 4 \log{\left (3 \right )}\\- 6 \log{\left (3 \right )} + 24\end{matrix}\right]$$
linearapprox(E**x,0,1)
$$\left[\begin{matrix}-10 + 4 e\\- 6 e + 18\end{matrix}\right]$$
from sympy import pi
tk=x*pi
linearapprox(cos(tk),0,1)
$$\left[\begin{matrix}\frac{12}{\pi^{2}}\\- \frac{24}{\pi^{2}}\end{matrix}\right]$$
linearapprox(log(x),1,2)
$$\left[\begin{matrix}- 4 \log{\left (2 \right )} + \frac{1}{2}\\-3 + 12 \log{\left (2 \right )}\end{matrix}\right]$$
Q17
Known experiment data is listed below:
| $x_i$ | 19 | 25 | 31 | 38 | 44 |
|---|---|---|---|---|---|
| $y_i$ | 19.0 | 32.3 | 49.0 | 73.3 | 97.8 |
Calculate a formula in the form of $y=a+bx^2$ for the data using least square method and its mean square error.
from scipy.optimize import curve_fit
import numpy as np
def dat(p,a,b):
return a+b*np.square(p)
fitpars,covmat = curve_fit(dat,[19,25,31,38,44],[19.0,32.3,49.0,73.3,97.8])
print fitpars, "MSE=", covmat.trace()
[ 0.97257866 0.05003512] MSE= 0.00454902193875
thus the solution is $y=0.973+0.05x^2$.
import matplotlib.pyplot as plt
fig, ax1 = plt.subplots()
f0=0.973 + 0.05*x**2
span= linspace(19,44,200)
ax1.plot(span, lambdify(x,f0)(span), lw=2, color="blue")
ax1.scatter([19,25,31,38,44],[19.0,32.3,49.0,73.3,97.8], lw=2, color="blue")
ax1.set_ylabel(r"Value", fontsize=18, color="blue")
for label in ax1.get_yticklabels():
label.set_color("blue")
PQ1
Evaluate the $f(x)=\frac{1}{1+25x^2}$ in $[-1,1]$ on points $x_i=-1+0.2i,i\in Z$. Calculate the cubic curve fit of the points above, and compare with the results in Chapter 2.
x0,y=[],[]
f0=1/(1+25*x**2)
for i in range(0,11):
x0.append(-1+0.2*i)
y.append( f0.subs(x,-1+0.2*i))
x0
$$\left [ -1.0, \quad -0.8, \quad -0.6, \quad -0.4, \quad -0.2, \quad 0.0, \quad 0.2, \quad 0.4, \quad 0.6, \quad 0.8, \quad 1.0\right ]$$
import sympy as syp
a,b,c,d=0,0,0,0
def data2(p,a,b,c,d):
return a+b*np.power(p,1)+c*np.power(p,2)+d*np.power(p,3)
coef_fit,covmat = curve_fit(data2,[-1.0,-0.8,-0.6,-0.4,-0.2,0.0,0.2,0.4,0.6,0.8,1.0],[0.0384615384615385,0.0588235294117647,0.1,0.2,0.5,1,0.5,0.2,0.1,0.0588235294117647,0.0384615384615385])
print coef_fit, "MSE=", covmat.trace()
[ 4.84124918e-01 -2.20984575e-08 -5.75182739e-01 1.57732056e-08] MSE= 0.295257692053
f_1=coef_fit[0]+coef_fit[1]*x+coef_fit[2]*x**2+coef_fit[3]*x**3
f_1
$$1.57732056305093 \cdot 10^{-8} x^{3} - 0.575182738803699 x^{2} - 2.20984575083349 \cdot 10^{-8} x + 0.484124918275589$$
plot(f0,f_1,(x,-1,1))
<sympy.plotting.plot.Plot at 0x7f74dcaf14d0>
In comparision with polynomial interpolation:
x1=linspace(-1,1,10)
y1=[]
f0=1/(1+25*x**2)
for i in x1:
y1.append( f0.subs(x,i))
y1
$$\left [ 0.0384615384615385, \quad 0.0620214395099541, \quad 0.114730878186969, \quad 0.264705882352941, \quad 0.764150943396226, \quad 0.764150943396226, \quad 0.264705882352941, \quad 0.114730878186969, \quad 0.0620214395099541, \quad 0.0384615384615385\right ]$$
x2=linspace(-1,1,20)
y2=[]
f0=1/(1+25*x**2)
for i in x2:
y2.append( f0.subs(x,i))
y2
$$\left [ 0.0384615384615385, \quad 0.0475876614816768, \quad 0.0603073838957568, \quad 0.0787178368948975, \quad 0.106615475487301, \quad 0.151299245599329, \quad 0.227616645649432, \quad 0.366125760649087, \quad 0.616040955631399, \quad 0.935233160621761, \quad 0.935233160621762, \quad 0.616040955631399, \quad 0.366125760649087, \quad 0.227616645649433, \quad 0.151299245599329, \quad 0.106615475487301, \quad 0.0787178368948975, \quad 0.0603073838957568, \quad 0.0475876614816768, \quad 0.0384615384615385\right ]$$
An routine to calculate the best interpolating quadratic
Now we will try to make a routine to interpolate a function using SymPy.
def make_standard_base(n):
x=Symbol('x')
bases=[]
matrix=[]
for i in range(0,n):
bases.append(x**i)
return bases
def make_equations(x0,x1,n):
x=Symbol('x')
bases=[]
matrix=[]
for i in range(0,n):
bases.append(x**i)
for i in range(0,n):
line=[]
for j in range(0,n):
line.append(integrate(bases[i]*bases[j],(x,x0,x1)))
matrix.append(line)
return Matrix(matrix)
def interpolate_with_base(equ,base, f,n,x0,x1):
cof=[]
for i in range(0,n):
cof.append(integrate(base[i]*f,(x,x0,x1)))
return equ.LUsolve(Matrix(n,1,cof))
make_equations(0,1,3)
$$\left[\begin{matrix}1 & \frac{1}{2} & \frac{1}{3}\\\frac{1}{2} & \frac{1}{3} & \frac{1}{4}\\\frac{1}{3} & \frac{1}{4} & \frac{1}{5}\end{matrix}\right]$$
interpolate_with_base(make_equations(0,1,3),make_standard_base(3),E**x,3,0,1)
$$\left[\begin{matrix}-105 + 39 e\\- 216 e + 588\\-570 + 210 e\end{matrix}\right]$$