数值分析方法 课件 3-3 Python程序在数值逼近中的应用

数值分析方法主编

李冬果李林高磊面向“四新”人才培养普通高等教育系列教材第三章数值逼近基础3.1插值逼近3.2曲线拟合3.3Python程序在数值逼近中的应用

目录/Contents

3.3Python程序在数值逼近中的应用差值算法Python实验拉格朗日（Lagrange）插值法Newton差值法的实现拟合算法Python实验线性拟合多项式拟合Python中的插值函数Numpy工具包中interp、ployfit、polyval等SciPy工具包中interpolate、curve_fit等拉格朗日（Lagrange）插值法deflagrange_interp(x_data,y_data,x_value):

n=len(x_data)

y=0

foriinrange(n):

p=1

forjinrange(i):

p=p*(x_value-x_data[j])/(x_data[i]-x_data[j])

forjinrange(i+1,n):

p=p*(x_value-x_data[j])/(x_data[i]-x_data[j])

y=y+y_data[i]*p

returny参数x_data和y_data是插值点的自变量和函数值，x_value是待求差值函数值的自变量值，为了同步计算多个点的差值函数值，x_value可以使用NumPy数组importnumpyasnpx_data=np.array([100,121,144],dtype=np.float64)y_data=np.sqrt(x_data)lagrange_interp(x_data,y_data,np.array([115]))结果为：array([10.72275551])对比真值：>>>np.sqrt(115)10.723805294763608defget_divided_differences(x,y):

n=len(x)

dd_matrix=np.zeros((n,n),dtype=float)

dd_matrix[:,0]=y

forjinrange(1,n):

foriinrange(j,n):

dd_matrix[i,j]=\

(dd_matrix[i,j-1]-dd_matrix[i-1,j-1])/(x[i]-x[i-j])

returndd_matrixNewton插值法函数中的符号“\”是Python中的换行提示符首先编写差商表计算函数>>>f=lambdax:x**7+x**4+3*x+1>>>x_data=[2**iforiinrange(9)]>>>x_data=np.array(x_data,dtype=np.float64)>>>y_data=f(x_data)>>>dd_matrix=get_divided_differences(x_data,y_data)>>>dd_matrix[7,7]1.0>>>dd_matrix[8,8]0.0编写Newton差值法的函数Newton插值法defnewton_interp(x_data,y_data,x_value):

n=len(x_data)

dd_matrix=get_divided_differences(x_data,y_data)

p=dd_matrix[n-1,n-1]

forkinrange(n-2,-1,-1):

p=dd_matrix[k,k]+(x_value-x_data[k])*p

returnp>>>x_data=np.array([100,121,144,169],dtype=np.float64)>>>y_data=np.sqrt(x_data)>>>newton_interp(x_data,y_data,115)10.723574251835121Numpy工具包中有一个基本的差值函数interp，而SciPy工具包更是包含一个专门的子包interpolate用于插值计算。这些包含了多种分段插值方法的实现，如常用的分段线性插值函数，三次样条插值类等。插值函数的使用

解：本题可以由如下代码完成：importnumpyasnpfromscipyimportinterpolatex=np.linspace(0,2*np.pi,9)y=np.sin(x)spl=interpolate.CubicSpline(x,y)x_value=np.linspace(0,2*np.pi,100)y_value_spl=spl(x_value)y_value_lin=erp(x_value,x,y)spl是由erpolate.CubicSpline和节点数据计算出的差值函数y_value_spl即为三次样条插值的结果y_value_lin则是由函数numpy.interp计算出的线性插值结果。通过matplotlib绘出插值图像：importmatplotlib.pyplotaspltplt.figure(figsize=(8,6),dpi=600)plt.plot(x,y,'kd')plt.plot(x_value,y_value_lin,'k:')plt.plot(x_value,y_value_spl,'k-')plt.legend(['Nodes','Linear','CubicSpline'])plt.show()线性拟合deflin_fit(x,y):

n=len(x)

s_x=sum(x)

s_y=sum(y)

s_xy=sum(x*y)

s_x2=sum(x**2)

a0=(s_x2*s_y-s_x*s_xy)/(n*s_x2-s_x**2)

a1=(n*s_xy-s_x*s_y)/(n*s_x2-s_x**2)

returna1,a0本例中涉及NumPy数组的运算，参数中的x和y须是数组才能正确运行例

试根据以下实验数据做拟合函数>>>x_data=np.array([1.36,1.37,1.95,2.28])>>>y_data=np.array([14.094,16.844,18.475,20.963])>>>lin_fit(x_data,y_data)(5.853760129659651,7.408457374392223)多项式拟合importnumpyasnpfromscipyimportlinalgasla

defpoly_fit(x,y,n):

n=n+1

s_x_list=np.array([sum(x**i)foriinrange(2*n)])

s_xy_list=np.array([sum(x**i*y)foriinrange(n)])

A=np.zeros((n,n))

foriinrange(n):

A[i,:]=s_x_list[i:(i+n)]

s_xy_list.reshape((n,1))

a_list=la.inv(A)@s_xy_list

returna_list[::-1]为了简便计算，在函数中直接调用的scipy.linalg.inv函数，并用NumPy矩阵运算计算了线性方程组的解。最后返回结果按照降幂排列例

试根据以下实验数据做拟合函数>>>x_data=np.array([0,1,2,3,4,5],dtype=np.float64)>>>y_data=np.array([5,2,1,1,2,3],dtype=np.float64)>>>poly_fit(x_data,y_data,2)array([0.5

,-2.78571429,

4.71428571])拟合函数的使用在NumPy工具包中提供了polyfit函数，用于一般的多项式拟合；NumPy还提供了函数polyval用于求系数向量形式的多项式函数值例

试根据以下实验数据做拟合函数>>>np.polyfit(x_data,y_data,2)array([0.5

,-2.78571429,

4.71428571])在SciPy工具包的优化子包optimize中，还包含有curve_fit函数，可求解一般曲线拟合问题拟合函数的使用fromscipy.optimizeimportcurve_fit

deffunc(t,a,k):

returna*np.exp(-k*t)

t_data=np.array([1,2,3,4,6,8,10],dtype=np.float64)c_data=np.array([2.0405,1.9050,1.7695,1.6339,1.3627,1.0917,0.8202],dtype=np.float64)popt,pcov=curve_fit(func,t_data,c_data)>>>print(popt)[2.301740630.09349889]即拟合后模型为：绘图以查看拟合效果：