- 我们在分析一般的复杂连续函数时,我们或许没办法直接做公式推导与数值计算,那么我们会寻求一个多项式来逼近这个函数,来简化我们的计算。常见的有连续函数的泰勒展开,今天我们来学习做数值计算时经常用到的——连续函数的最佳平方逼近。
定义1.内积与内积空间
- 在高等代数中我们曾经学习到了向量空间的内积,在这里我们将会定义函数的内积与内积空间。
设U是实线性空间,是定义在区间上的连续函数,在U上定义一个二元实值函数
其中为上的一个权函数,则称为定义在空间U上的内积,若满足以下性质:
- ,即满足可交换性。
- ,,即满足齐次性。
- ,即满足可加性。
- ,,当且仅当,即满足非负性。
则称U为[a.b]上的一个内积空间。
定义2.内积赋范空间
- 设U是一个内积空间,有定义
为U上的一个范数,称U为一个赋范空间,且内积范数满足:
定义3.最佳平方逼近
- 在我们学习了内积和内积空间的基础知识之后,我们要利用内积这一工具来开展接下来的工作,一般来讲,我们会选用距离最近来描述一个拟合的效果好坏程度,在内积空间中,我们通常会使用2-范数来描述距离。
对于函数且中的一组子集,若存在使得
则称是在子集 中的最佳平方逼近函数。
求解
- 从上边我们可以知道,若为了使得最小,则一般地,我们会采用求偏导数为零的点,来计算最值,即,通过化简我们可以得到如下的方程:
那么可以得到一个线性方程组:
以上的方程组称为正规方程组。
若取,则若要求n次最佳平方逼近多项式
此时,
用H表示对于的矩阵,即
称H为希尔伯特矩阵,记则
的解即为所求。
勒让德多项式
- 当区间为[0,1],权函数时,由正交化得到的多项式成为勒让德多项式,并用表示,罗德利克给出了勒让德多项式德简单表达式:
用python求解
- 我们要求已知的函数在上的的最佳平方逼近,我们按照本节课所学的知识给出我们
- 导入必要的库
import numpy as np
import math
import sympy
from scipy import integrate
import matplotlib.pyplot as plt
plt.style.use("ggplot")
- 定义求解的函数
def function_d(a,b,n):#求d_k,k=0,1,...,n
x = sympy.symbols('x')
f = sympy.sin(x)
G = [sympy.Symbol('t') for c in range(n)]
v = np.zeros(n)
for i in range(n):
G[i]=f*pow(x,i)
v[i]=sympy.integrate(G[i],(x,a,b))
return v
- 定义求解正规矩阵H函数
def function_H(a,b,n):#求解正规矩阵H
y = np.zeros((n,n))
for i in range(n):
for j in range(n):
y[i,j]=1
y[i][j]=integrate.quad(lambda x:pow(x,i+j), a, b)[0]
return y
- 定义求解系数向量a的函数
def function_a(a,b,n):
q = np.zeros((n,n))
g = np.zeros((n, n))
q=function_H(a,b,n)#求解正规方程组H
g=np.linalg.inv(q)
s=v=np.zeros(n)
s=function_d(a,b,n)#求d_k
v=np.dot(g,s)#a=H^-1*d
return v
- 输入逼近函数的次数并求解表达式,这里我们选用三次多项式进行函数逼近。
m=int(input("请输入次数:"))
z=function_a(0,math.pi,m+1)#求a:下限,b:下限,n:n-1次多项式
x = sympy.symbols('x')
s_0 = z[0]
for i in range(m+1):
if i<m:
s_0 =s_0+pow(x,i+1)*z[i+1]
else:
break
print ("s_0={}".format(s_0))
这里求解的到表达式
- 绘制和逼近函数的图像
def s_1(x):
s_x = z[0]
for i in range(m+1):
if i<m:
s_x =s_x+pow(x,i+1)*z[i+1]
else:
break
return s_x
x = np.linspace(a,b,10000)
y = np.sin(x)
plt.style.use("ggplot")
plt.plot(x, s_1(x), ls="-",color="r",marker=",",lw=2, label='s_0(x)')
plt.plot(x, y, ls="-",color="b",marker=",",lw=2, label='y=sin(x)')
plt.xlabel("X")
plt.ylabel("y")
plt.title("figure")
plt.legend()
plt.show()
得到如下图像
这里可以发现,我们所做的逼近函数和相当接近。
参考资料:
《数值分析》第5版,清华大学出版社
参考链接: