# 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:

1. $f(x)=\frac{1}{x},[1,3]$
2. $f(x)=e^x,[0,1]$
3. $f(x)=\text{cos}\pi x,[0,1]$
4. $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]$$