连续函数的最佳平方逼近

最新推荐文章于 2024-01-09 02:03:15 发布

m0_61779914

最新推荐文章于 2024-01-09 02:03:15 发布

阅读量3.8k

点赞数 35

本文链接：https://blog.csdn.net/m0_61779914/article/details/125341803

版权

最佳平方逼近内积空间勒让德多项式数值计算 Python实现

关键词由CSDN通过智能技术生成

我们在分析一般的复杂连续函数时，我们或许没办法直接做公式推导与数值计算，那么我们会寻求一个多项式来逼近这个函数，来简化我们的计算。常见的有连续函数的泰勒展开，今天我们来学习做数值计算时经常用到的——连续函数的最佳平方逼近。

定义1.内积与内积空间

在高等代数中我们曾经学习到了向量空间的内积，在这里我们将会定义函数的内积与内积空间。

设U是实线性空间， $f(x),g(x)$ 是定义在区间 $\left [ a ,b\right ]$ 上的连续函数,在U上定义一个二元实值函数

$< f(x),g(x)> = \int_{a}^{b}w(x)f(x)g(x)dx$

其中 $w(x)$ 为 $\left [ a ,b\right ]$ 上的一个权函数，则称为 $f(x),g(x)$ 定义在空间U上的内积，若满足以下性质：

$<f,g>=<g,f>,\forall f,g\in U$ ，即满足可交换性。
$<kf,g>=k<g,f>,\forall f,g\in U$ ， $\forall k\in \mathbb{R}$ ，即满足齐次性。
$<f+g,g>=<f,g>+<g,g>,\forall f,g\in U$ ,即满足可加性。
$<f,f>\geq0,\forall f \in \mathbb{R}$ ， $<f,f>=0$ ，当且仅当 $f\equiv 0$ ，即满足非负性。

则称U为[a.b]上的一个内积空间。

定义2.内积赋范空间

设U是一个内积空间，有定义

$\left \| f \right \|=\sqrt{<f,f>},\forall f \in U$

为U上的一个范数，称U为一个赋范空间，且内积范数满足:

$\left \| f \pm g \right \|^{2}=\left \| f \right \|^{2} \pm 2<f,g> +\left \| g \right \|^{2},\forall f,g \in U$
$\left | <f,g> \right |\leq \left \| f\left \| \right \| g\right \|,\forall f,g\in U$
$\left \| f+g \right \|^{2}+\left \| f-g \right \|^{2}=2\left \| f \right \|^{2}+2\left \| g \right \|^{2},\forall f,g \in U$

定义3.最佳平方逼近

在我们学习了内积和内积空间的基础知识之后，我们要利用内积这一工具来开展接下来的工作，一般来讲，我们会选用距离最近来描述一个拟合的效果好坏程度,在内积空间中，我们通常会使用2-范数来描述距离。

对于函数 $f= f(x)\in \mathbb{C}\left[a,b \right ]$ 且 $\mathbb{C}\left [a,b \right ]$ 中的一组子集 $\varphi =span\left \{ \varphi_{0}(x),\varphi_{1}(x),...,\varphi_{n}(x) \right \}$ ,若存在 $S^{\ast }(x)\in \varphi$ 使得

$\left \| f(x) -S_{0 }(x)\right \|_{2}^{2}=\min_{S_{0}(x)\in \varphi }\left \| f(x) -S_{0 }(x)\right \|_{2}^{2}$

$=\min_{S_{0}(x)\in \varphi }\int_{a}^{b}w(x)[f(x)-S_{0}(x)]^{2}dx$

则称 $S_{0}(x)$ 是 $f= f(x)\in \mathbb{C}\left[a,b \right ]$ 在子集 $\varphi \in \mathbb{C}\left [a,b \right ]$ 中的最佳平方逼近函数。

求解

从上边我们可以知道 $S_{0}(x)=a_{0}\varphi _{0}(x)+a_{1}\varphi _{1}(x)+...+a_{n}\varphi _{n}(x)$ ，若为了使得 $\left \| f(x) -S_{0 }(x)\right \|_{2}^{2}$ 最小，则一般地，我们会采用求偏导数为零的点，来计算最值，即 $\frac{\partial S_{0}}{\partial a_{j}}=0,j\in 0,1,...,n$ ,通过化简我们可以得到如下的方程:

$\sum_{i=1}^{n}a_{j}\int_{a}^{b}w(x)\varphi _{i}(x)\varphi _{j}(x)dx=\int_{a}^{b}w(x)f(x)\varphi _{j}(x)dx$

那么可以得到一个线性方程组：

$\begin{bmatrix} <\varphi _{0}(x), \varphi _{0}(x)> & <\varphi _{0}(x), \varphi _{1}(x)> & ... & <\varphi _{0}(x), \varphi _{n}(x)>\\ <\varphi _{1}(x), \varphi _{0}(x)>& <\varphi _{1}(x), \varphi _{1}(x)> & ... & <\varphi _{1}(x), \varphi _{n}(x)> \\ \vdots & \vdots & &\vdots \\ <\varphi _{n}(x), \varphi _{0}(x)>& <\varphi _{n}(x), \varphi _{1}(x)> & ... & <\varphi _{n}(x), \varphi _{n}(x)> \end{bmatrix}\begin{pmatrix} a_{0}\\ a_{1}\\ \vdots \\ a_{n} \end{pmatrix}=\begin{pmatrix} < f(x),\varphi _{0}(x)> \\ < f(x),\varphi _{1}(x)> \\ \vdots \\ < f(x),\varphi _{n}(x)> \end{pmatrix}$

以上的方程组称为正规方程组。

若取 $\varphi _{i}(x)=x^{k},w (x)\equiv 1,f(x)\in \mathbb{C}[0,1]$ ,则若要求n次最佳平方逼近多项式

$S_{0}(x)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n}x^{n}$

此时，

$<\varphi_{i}(x) ,\varphi_{j}(x)> =\int_{0}^{1}x^{i+j}dx=\frac{1 }{ i+j+1}$

$<f(x),\varphi_{i}(x)>=\int_{0}^{1}f(x)x^{k}dx\equiv d_{k}$

用H表示 $G_{n}=G(1,x,x^{2},...,x^{n})$ 对于的矩阵，即

$H=\begin{pmatrix} 1 & \frac{1}{2}& ... & \frac{1}{n+1}\\ \frac{1}{2} & \frac{1}{3} & ... &\frac{1}{n+2} \\ \vdots & \vdots & &\vdots\\ \frac{1}{n+1} & \frac{1}{n+2} & ... & \frac{1}{2n} \end{pmatrix}$

称H为希尔伯特矩阵，记 $a=(a_{0},a_{1},...,a_{n})^{T},d=(d_{0},d_{1},...,d_{n})^{T},$ 则

$Ha=d$

的解 $a_{i}=a_{i}^{*},(i=0,1,...,n)$ 即为所求。

勒让德多项式

当区间为[0,1],权函数 $w(x)\equiv 1$ 时，由 $\left \{1,x,...x^{n},... \right \}$ 正交化得到的多项式成为勒让德多项式，并用 $P_{i}(x)$ 表示，罗德利克给出了勒让德多项式德简单表达式:

$P_{x}=1,P_{n}(x)=\frac{1}{2n!}\frac{\mathrm{d^{n}} }{\mathrm{d} x^{n}}(x^{2}-1)^{n},n=1,2,...,n,...$

用python求解

我们要求已知的函数 $y= sin(x)$ 在 $[0,\pi]$ 上的的最佳平方逼近，我们按照本节课所学的知识给出我们
导入必要的库

import numpy as np
import math
import sympy
from scipy import integrate
import matplotlib.pyplot as plt
plt.style.use("ggplot")

定义求解 $d_{k}$ 的函数

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))

这里求解的到表达式 $s_{0}(x)=-1.066\times10^{-14}x^{3} - 0.418x^{2} + 1.312x - 0.051$

绘制 $y=sin(x)$ 和逼近函数 $s_{0}(x)=-1.066\times10^{-14}x^{3} - 0.418x^{2} + 1.312x - 0.051$ 的图像

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()

得到如下图像