数值分析实习作业（各种插值函数与积分公式的python代码实现）

最新推荐文章于 2024-04-30 19:24:04 发布

Angra-Mainyu

最新推荐文章于 2024-04-30 19:24:04 发布

阅读量1.2w

点赞数 26

文章标签： python 插值函数数值分析

本文链接：https://blog.csdn.net/qq_15800759/article/details/80629629

版权

数值分析实习作业

1.对函数 $f(x)$ 进行插值

f (x) = 1 1 + x 2

$f(x)=\frac{1}{1+x^{2}}$

(1)令插值节点为等距节点 $-5,-4,-3,-2,-1,0,1,2,3,4,5$ ，在这些节点处对 $f(x)$ 进行Lagrange插值和Hermite插值；

Lagrange 插值：

算法描述：

对某个多项式函数，已知有给定的k + 1个取值点：

(x 0, y 0), \dots, (x k, y k)

$(x_{0},y_{0}),\ldots ,(x_{k},y_{k})$
其中

xj x j $x_{j}$ 对应着自变量的位置，而

yj y j $y_{j}$ 对应着函数在这个位置的取值。

假设任意两个不同的xj都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：

L (x) = \sum j = 0 k y j l j (x)

$L(x)=\sum_{j=0}^ky_{j}l_{j}(x)$
其中每个

li(x) l i ( x ) $l_i(x)$ 被称为拉格朗日基本多项式，也叫 拉格朗日插值基函数：

l j (x) = \prod i = 0, i \neq j k x - x i x j - x i

$l_{j}(x)=\prod_{i=0,i\neq{j}}^k \frac{x-x_{i}}{x_{j}-x_{i}}$

效果图：

这里写图片描述
可见在函数两端有明显的龙格现象，验证了高次多项式插值函数存在的这个问题

在数值分析领域中，龙格现象是在一组等间插值点上使用具有高次多项式的多项式插值时出现的区间边缘处的振荡问题。它是由卡尔·龙格（Runge）在探索使用多项式插值逼近某些函数时的错误行为时发现的。它表明使用高次多项式插值并不总能提高准确性。该现象与傅里叶级数近似中的吉布斯现象相似。

平均误差值:

计算区间 $[-5,5]$ 内 $101$ 个等距分布点处的误差绝对值的平均值

0.285755407117

实现代码如下：

# 生成Lagrange插值函数，返回以xi做参数的函数指针
def Lagrange(x, f):
    y = f(x)
    def lx(xi):
        result = 0.0
        for i in range(len(x)):
            tmp = y[i]
            for j in range(len(y)):
                if i != j:
                    tmp *= (xi - x[j]) / (x[i] - x[j])
            result += tmp
        return result
    return lx

Hermite 插值：

算法描述：

埃尔米特插值是另一类插值问题，这类插值在给定的节点处，不但要求插值多项式的函数值与原函数值相同。同时还要求在节点处，插值多项式的一阶直至指定阶的导数值，也与被插函数的相应阶导数值相等，这样的插值称为埃尔米特(Hermite)插值。

Hermite插值在不同的节点，提出的插值条件个数可以不同，若在某节点 $x_{i}$ ，要求插值函数多项式的函数值，一阶导数值，直至 $m_{i}-1$ 阶导数值均与被插函数的函数值相同及相应的导数值相等。我们称 $x_{i}$ 为 $m_{i}$ 重插值点节,因此，Hermite插值应给出两组数，一组为插值点 $\{x_{i}\}_{i=0}^{n}$ 节点，另一组为相应的重数标号 $\{m_{i}\}_{i=0}^{n}$ 。
若 $\sum _{i=0}^{n}m_{i}=N+1$ ，这就说明了给出的插值条件有 $N+1$ 个，为了保证插值多项式的存在唯一性，这时的Hermite插值多项式应在 $P_{n}$ 上求得.

当每一个插值节点再有一个对于导数的约束条件时，此时线性方程变成了2*（n + 1）= 2n + 2个，也就是需要构造一个不超过2n+1次的多项式。

构造出的埃尔米特插值函数应该是这样一个形式：

H 2 n + 1 (x) = \sum i = 0 n f (x i) α i (x) + \sum i = 0 n f' (x i) β i (x)

$H_{2n+1}(x) =\sum_{i=0}^{n}f(x_{i})\alpha _{i}(x)+\sum_{i=0}^{n}f{}'(x_{i})\beta _{i}(x)$

两个基函数：

{α i (x j) = δ i j, α' i (x j) = 0, j = 0, 1, 2, . . ., n, {β i (x j) = 0, β' i (x j) = δ i j, j = 0, 1, 2, . . ., n, δ i j = {1, j = i, 0, j \neq i, j = 0, 1, 2, . . ., n .

$\begin{array}{l}\left\{\begin{array}{l}\alpha_i(x_j)=\delta_{ij},\\\alpha'_i(x_j)=0,\;\end{array}\right.j=0,1,2,...,n,\\\\\left\{\begin{array}{l}\beta_i(x_j)=0,\\\beta'_i(x_j)=\delta_{ij},\end{array}\right.j=0,1,2,...,n,\\\\\delta_{ij}=\left\{\begin{array}{l}1,j=i,\\0,j\neq i,\end{array}j=0,1,2,...,n.\right.\end{array}$

根据课本上提供的公式构造出两个符合条件的基函数 $\alpha_i\left(x\right)$ 和 $\beta_i\left(x\right)$ 如下：

α i (x) = l 2 i (x) [1 + 2 (x i - x) \sum k = 0, k \neq i n 1 x i - x k]

$\alpha_i(x)=l_i^2(x)\;\lbrack1+2(x_i-x)\sum_{k=0,k\neq i}^n\frac1{x_i-x_k}\rbrack$

β i (x) = (x - x i) l 2 i (x)

$\beta_i(x)=(x-x_i)l_i^2(x)$

其中 $l_i(x)$ 是拉格朗日插值基函数

l j (x) = \prod i = 0, i \neq j k x - x i x j - x i

$l_{j}(x)=\prod_{i=0,i\neq{j}}^k \frac{x-x_{i}}{x_{j}-x_{i}}$
因此我们只需要根据公式计算出

αi(x) α i ( x ) $\alpha_i(x)$ 和

βi(x) β i ( x ) $\beta_i(x)$ 这两个基函数，就可以很简单地得出赫尔米特插值函数了

效果图如下：

这里写图片描述
可见此时Hermite插值函数依旧有着很严重的龙格现象

平均误差值:

计算区间[-5,5]内 101 个等距分布点处的误差绝对值的平均值：

0.350297537083

实现代码如下：

根据公式 $l_{j}(x)=\prod_{i=0,i\neq{j}}^k \frac{x-x_{i}}{x_{j}-x_{i}}$ 生成朗格朗日基函数:

# 生成朗格朗日基函数
def get_li(i, x_all):
    def li(xx):
        a = 1.0
        b = 1.0
        for j in range(len(x_all)):
            if j == i:
                continue
            a *= (xx-x_all[j])
            b *= (x_all[i]-x_all[j])
        return a / b
    return li

根据公式 $\alpha_i(x)=l_i^2(x)\;\lbrack1+2(x_i-x)\sum_{k=0,k\neq i}^n\frac1{x_i-x_k}\rbrack$ 生成 $\alpha_i(x)$ 基函数：

# 求alpha基函数
def get_alpha(i, x_all=[]):
    def alpha(xx):
        result = 0.0
        for k in range(len(x_all)):
            if k == i:
                continue
            result += 1.0/(x_all[i]-x_all[k])
        result = (1 + 2*(x_all[i]-xx) * result)
        li = get_li(i, x_all)
        result *= li(xx)**2
        return result
    return alpha

根据公式 $\beta_i(x)=(x-x_i)l_i^2(x)$ 生成 $\beta_i(x)$ 基函数：

# 求alpha基函数
def get_beta(i, x_all=[]):
    def beta(xx):
        li = get_li(i, x_all)
        result = (xx - x_all[i]) * li(xx)**2
        return result
    return beta

代入 $H_{2n+1}(x) =\sum_{i=0}^{n}f(x_{i})\alpha _{i}(x)+\sum_{i=0}^{n}f{}'(x_{i})\beta _{i}(x)$ 生成Hermite插值函数：

def Hermite(x, f, f_):
    def He(xi):
        result = 0.0
        n = len(x)
        for i in range(n):
            ai = get_alpha(i, x)
            bi = get_beta(i, x)
            result += y[i] * ai(xi) + y_[i] * bi(xi)
        return result
    return He

(2)令插值节点为区间 $[-5,5]$ 上的 $11$ 次切比雪夫多项式的零点，在这些节点处对f 进行 Lagrange 插值；

算法描述：

其实该题的意思就是使用切比雪夫多项式零点插值来处理这个函数。
第一类切比雪夫多项式由以下三角恒等式确定:

T n (cos (θ)) = cos (n θ)

$T_{n}(\cos(\theta ))=\cos(n\theta )$
其中

n=0,1,2,3,.....cosnθ n = 0 , 1 , 2 , 3 , . . . . . cos ⁡ n θ $n = 0, 1, 2, 3, .... . \cos n\theta$ , 是关于

cosθ cos ⁡ θ $\cos\theta$ 的

n n $n$ 次多项式,
用显式来表示：

T_{n} (x) = {\begin{cases} \cos (n \arccos (x)), & x \in [- 1, 1] \\ \cosh (n a r c c o s h (x)), & x \geq 1 \\ (- 1)^{n} \cosh (n a r c c o s h (- x)), & x \leq - 1 \end{cases}

$T_{n}(x)={\begin{cases}\cos(n\arccos(x)),&\ x\in [-1,1]\\\cosh(n\,\mathrm {arccosh} (x)),&\ x\geq 1\\(-1)^{n}\cosh(n\,\mathrm {arccosh} (-x)),&\ x\leq -1\\\end{cases}}$

$n$ 次切比雪夫多项式在区间[−1,1]上都有 $n$ 个不同的根, 称为切比雪夫根, 有时亦称做切比雪夫节点，因为是多项式插值时的插值点 . 从三角形式中可看出 $T_n$ 的 $n$ 个根分别是：

x_{i} = \cos (\frac{2 i - 1}{2 n} π), i = 1, . . ., n

$x_i=\cos(\frac{2i-1}{2n}\pi),i=1,...,n$

根据题目条件计算出11个切比雪夫多项式的零点：
[ 4.949107, 4.548160, 3.778748, 2.703204, 1.408663, 0.000000, -1.408663, -2.703204, -3.778748, -4.548160, -4.949107]
接着以这些节点作为朗格朗日插值节点进行拉格朗日插值即可

效果图如下：

这里写图片描述
函数两端的波动相较之前两个插值函数已经小了很多

平均误差值:

0.0486976071887

实现代码如下：

根据公式 $x_i=\cos(\frac{2i-1}{22}\pi),i=1,...,11$ 获取11次切比雪夫多项式的零点

# 获取11次切比雪夫多项式的零点
def ChebyshevX():
    n = 11
    x = np.arange(n) + 1.0
    # print x
    xi = np.cos((2*x-1)/(2*n)*np.pi)
    return xi * 5

然后使用这些节点调用上一问已经完成的函数进行拉格朗日插值:

# 获取切比雪夫多项式零点来代替原插值节点
x = ChebyshevX()
# 获取切比雪夫零点插值函数
Lx = Lagrange(x, f)

得到的函数Lx即为切比雪夫多项式零点插值函数

(3)令插值节点为等距节点 $-5,-4,-3,-2,-1,0,1,2,3,4,5,$ 在这些节点处对 $f(x)$ 进行分段线性插值和分段三次 Hermite 插值；

分段线性插值：

算法描述：

分段线性插值就是通过插值点用折线段连接起来逼近 $f(x)$ .
设已知节点 $a=x_0<x_1<...<x_n= b$ 上的函数值 $f_0,f_1,...,f_n$ ，记 $h_k = x_{k+1}-x_k$ ， $h= max h_k$ ，求一折线函数 $I_h(x)$ 满足：
(1). $\begin{array}{l}I_h(x)\in C\lbrack a,b\rbrack\\\end{array}$
(2). $\begin{array}{l}I_h(x)=f_k(k=0,1,...,n)\\\end{array}$
(3). $\begin{array}{l}I_h\left(x\right)\\\end{array}$ 在每个小区间 $[x_k,x_{k+1}]$ 上是线性函数
则称 $I_h\left(x\right)$ 为分段线性插值函数
根据定义可知对于每个 $I_h(x)$ 在每个小区间 $[x_k,x_{k+1}]$ 上可以表示为：

I h (x) = x - x k + 1 x k - x k + 1 f k + x - x k x k + 1 - x k f k + 1

$I_h(x)\;=\frac{x-x_{k+1}}{x_k-x_{k+1}}f_k+\frac{x-x_k}{x_{k+1}-x_k}f_{k+1}$
其中：

x k \leq x \leq x k + 1 k = 0, 1, . . ., n - 1

$\begin{array}{l}x_k\leq x\leq x_{k+1}\\k=0,1,...,n-1\end{array}$
在本题中只需要将

x x $x$ 对应在每个区间

[x_{k}, x_{k + 1}]

$[x_k,x_{k+1}]$ 按照公式求出即可

效果图如下：

这里写图片描述
可见分段线性插值虽然插值的结果不光滑，但是在速度上和误差上都表现得比较好，不会出现龙格现象

平均误差值：

0.0147779905191

实现代码如下：

# 传入插值节点序列x和原函数f，获取分段线性插值函数
def subLiner(xi, f):
    yi = f(xi)

    def Lh(x):
        result = 0.0
        # 利用k来定位x所处的区间
        k = 0
        while k < len(xi) and not (xi[k] <= x and xi[k+1] >= x):
            k += 1
        result = (x - xi[k+1])/(xi[k] - xi[k+1]) * yi[k] + \
            (x - xi[k])/(xi[k+1] - xi[k]) * yi[k+1]
        return result
    return Lh

分段三次埃尔米特插值：

算法描述：

分段线性插值函数 $I_h(x)$ 的导数是间断的，若在节点 $x_k(k=0,1,...,n)$ 上除已知函数值 $f_k$ 之外还给出导数值 $f'_k=m_k(k=0,1,...,n)$ ，这样就可构造一个导数连续的分段插值函数 $I_h(x)$ ，称为分段三次埃尔米特插值函数，它满足条件：
(1). $\begin{array}{l}I_h(x)\in C'\lbrack a,b\rbrack\\\end{array}$
(2). $\begin{array}{l}I_h(x)=f_k,I'_h(x)=f'_k(k=0,1,...,n)\\\end{array}$
(3). $\begin{array}{l}I_h\left(x\right)\\\end{array}$ 在每个小区间 $[x_k,x_{k+1}]$ 上是三次多项式
根据两点三次插值多项式可知， $I_h(x)$ 在区间 $[x_k,x_{k+1}]$ 上的表达式为：

I h (x) = (x - x k + 1 x k - x k + 1) 2 (1 + 2 x - x k x k + 1 - x k) f k + (x - x k x k + 1 - x k) 2 (1 + 2 x - x k + 1 x k - x k + 1) f k + 1

$I_h\left(x\right)=\left(\frac{x-x_{k+1}}{x_k-x_{k+1}}\right)^2\left(1+2\frac{x-x_k}{{x_k}_{+1}-x_k}\right)f_k+\left(\frac{x-x_k}{x_{k+1}-x_k}\right)^2\left(1+2\frac{x-x_{k+1}}{{x_k}_{}-x_{k+1}}\right)f_{k+1}$

+ (x - x k + 1 x k - x k + 1) 2 (x - x k) f' k + (x - x k x k + 1 - x k) 2 (x - x k + 1) f' k + 1

$\;\;\;\;\;\;\;+\left(\frac{x-x_{k+1}}{x_k-x_{k+1}}\right)^2(x-x_k)f'_k+\left(\frac{x-x_k}{x_{k+1}-x_k}\right)^2(x-x_{k+1})f'_{k+1}$
上式对于

k=0,1,...,n−1 k = 0 , 1 , . . . , n − 1 $k=0,1,...,n-1$ 成立

效果图如下：

这里写图片描述
可见分段三次赫尔米特插值拟合性，光滑程度上都表现得非常优秀

平均误差值:

0.00132363164301

实现代码如下：

# 传入插值节点序列x、原函数f和导函数f_，生成分段埃尔米特插值函数
def subHermite(xi, f, f_):
    yi = f(xi)
    y_i = f_(xi)

    def lh(x):
        result = 0.0
        # 利用k来定位x所处的区间
        k = 0
        while k < len(xi) and not (xi[k] <= x and xi[k+1] >= x):
            k += 1
        result = ((x - xi[k+1]) / (xi[k] - xi[k+1]))**2 * (1 + 2* (x-xi[k])/ (xi[k+1]-xi[k])) * yi[k] +\
            ((x-xi[k])/(xi[k+1]-xi[k]))**2 * (1 + 2 * (x-xi[k+1])/(xi[k] - xi[k+1]) )*yi[k+1] + \
            ((x-xi[k+1])/(xi[k]-xi[k+1]))**2 * (x-xi[k]) * y_i[k] + \
            ((x-xi[k])/(xi[k+1]-xi[k]))**2 *(x-xi[k+1]) * y_i[k+1]
        return result
    return lh

(4)令插值节点为等距节点 $-5,-4,-3,-2,-1,0,1,2,3,4,5,$ 在这些节点处对 $f(x)$ 进行三次样条插值，其中边界条件为第一类边界条件，在两个端点处导数为 $1$ ，即 $f'(-5)=f'(5)=1$

第一类三次样条插值：

算法描述：

对于 $n+1$ 个给定点的数据集 ${x_i}$ ，我们可以用 $n$ 段三次多项式在数据点之间构建一个三次样条。如果

S (x) = {\begin{matrix} S_{0} (x), x \in [x_{0}, x_{1}] \\ S_{1} (x), x \in [x_{1}, x_{2}] \\ \dots \\ S_{n - 1} (x), x \in [x_{n - 1}, x_{n}] \end{matrix}

$S(x)=\left\{{\begin{matrix}S_{0}(x),\ x\in [x_{0},x_{1}]\\S_{1}(x),\ x\in [x_{1},x_{2}]\\\cdots \\S_{{n-1}}(x),\ x\in [x_{{n-1}},x_{n}]\end{matrix}}\right.$
表示对函数

f f $f$ 进行插值的样条函数，那么需要：

插值特性， $S(x_i)=f(x_i)$
样条相互连接， $S_{i-1}(x_i) = S_i(x_i), i=1,...,n-1$
两次连续可导， $S'_{i-1}(x_i) = S'_i(x_i)$ 以及 $S''_{i-1}(x_i) = S''_i(x_i), i=1,...,n-1.$

由于每个三次多项式需要四个条件才能确定曲线形状，所以对于组成S的 $n$ 个三次多项式来说，这就意味着需要 $4n$ 个条件才能确定这些多项式。但是，插值特性只给出了 $n + 1$ 个条件，内部数据点给出 $n+1-2=n-1$ 个条件，总计是 $4n-2$ 个条件。我们还需要另外两个条件，根据不同的因素我们可以使用不同的条件。

其中一项选择条件可以得到给定u与v的钳位三次样条，
$S'(x_{0})=u\,\!$
$S'(x_{k})=v\,\!$
另外，我们可以设
$S''(x_{0})=S''(x_{n})=0\,\!.$
这样就得到自然三次样条。自然三次样条几乎等同于样条设备生成的曲线。

在这些所有的二次连续可导函数中，钳位与自然三次样条可以得到相对于待插值函数f的最小震荡。

如果选择另外一些条件，
$S(x_{0})=S(x_{n})\,\!$
$S'(x_{0})=S'(x_{n})\,\!$
$S''(x_{0})=S''(x_{n})\,\!$
可以得到周期性的三次样条。

下面我们利用 $S(x)$ 的二阶导数值 $S''(x_j)=M_j(j=0,1,...,n)$ 表达 $S(x)$ .由于 $S(x)$ 在区间 $[x_j,x_{j+1}]$ 上是三次多项式，故 $S''(x)$ 在 $[x_j,x_{j+1}]$ 上是线性函数，可表示为

S'' (x) = M j x j + 1 - x h j + M j + 1 x - x j h j

$S''(x)=M_j\frac{x_{j+1}-x}{h_j}+M_{j+1}\frac{x-x_j}{h_j}$
对

S′′(x) S ″ ( x ) $S''(x)$ 积分两次并且利用

S(xj)=yj S ( x j ) = y j $S(x_j)=y_j$ 及

S(xj+1)=yj+1 S ( x j + 1 ) = y j + 1 $S(x_{j+1})=y_{j+1}$ ，可求出定积分常数，于是得到三次样条表达式

S (x) = M j ( x j + 1 - x ) 2 6 h j + M j + 1 ( x - x j ) 2 6 h j + (y j - M j h 2 j 6) x j + 1 - x h j + (y j + 1 - M j + 1 h 2 j 6) x - x j h j, j = 0, 1, . . ., n - 1

$\begin{array}{l}S(x)=M_j\frac{(x_{j+1}-x)^2}{6h_j}+M_{j+1}\frac{(x-x_j)^2}{6h_j}+\left(y_j-\frac{M_jh_j^2}6\right)\frac{x_{j+1}-x}{h_j}+\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left(y_{j+1}-\frac{M_{j+1}h_j^2}6\right)\frac{x-x_j}{h_j},\;\;\;j=0,1,...,n-1\end{array}$
在这里

Mj(j=0,1,...n) M j ( j = 0 , 1 , . . . n ) $M_j(j=0,1,...n)$ 是未知的.
利用公式求得：

μ j M j - 1 + 2 M j + λ j M j + 1 = d j, j = 1, 2, . . ., n - 1

$\mu_jM_{j-1}+2M_j+\lambda_jM_{j+1}=d_j,\;\;j=1,2,...,n-1$
其中

μ j = h j - 1 h j - 1 + h j, λ j = h j h j - 1 + h j,

$\begin{array}{l}\mu_j=\frac{h_{j-1}}{h_{j-1}+h_j},\;\;\lambda_j=\frac{h_j}{h_{j-1}+h_j},\\\end{array}$

d j = 6 f [ x j , x j + 1 ] - f [ x j - 1 , x j ] h j - 1 + h j = 6 f [x j - 1, x j, x j + 1], j = 1, 2, . . . n - 1

$d_j=6\frac{f\lbrack x_j,x_{j+1}\rbrack-f\lbrack x_{j-1},x_j\rbrack}{h_{j-1}+h_j}=6f\lbrack x_{j-1},x_j,x_{j+1}\rbrack,\;j=1,2,...n-1$
令

λ0=1,d0=6h0(f[x0,x1]−f′0),μn=1,dn=6hn−1(f′n−f[xn−1,xn]) λ 0 = 1 , d 0 = 6 h 0 ( f [ x 0 , x 1 ] − f 0 ′ ) , μ n = 1 , d n = 6 h n − 1 ( f n ′ − f [ x n − 1 , x n ] ) $\lambda_0=1,d_0=\frac6{h_0}(f\lbrack x0,x1\rbrack-f'_0),\;\mu_n=1,d_n=\frac6{h_n-1}(f'_n-f\lbrack x_{n-1},x_n\rbrack)$ ，然后可以列矩阵求解：

⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ 2 μ 2 λ n λ 1 2 ⋱ λ 2 ⋱ μ n - 1 ⋱ 2 μ n μ 1 λ n - 1 2 ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ M 1 M 2 ⋮ M n - 1 M n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ = ⎛ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ d 1 d 2 ⋮ d n - 1 d n ⎞ ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟

$\begin{pmatrix}2&\lambda_1&&&\mu_1\\\mu_2&2&\lambda_2&&\\&\ddots&\ddots&\ddots&\\&&\mu_{n-1}&2&\lambda_{n-1}\\\lambda_n&&&\mu_n&2\end{pmatrix}\begin{pmatrix}M_1\\M_2\\\vdots\\M_{n-1}\\M_n\end{pmatrix}=\begin{pmatrix}d_1\\d_2\\\vdots\\d_{n-1}\\d_n\end{pmatrix}$
解出

M M $M$ 之后代入公式即可解出

S (x)

$S(x)$

条件：

f' (- 5) = f' (5) = 1

$f'(-5)=f'(5)=1$

效果图如下：

这里写图片描述
因为在两个端点处要保证导数值为1，所以插值图像中两端点附近出现了一定程度的抖动

平均误差值:

0.0313803197288

实现代码如下：

# 传入插值节点序列x、原函数f、导函数f_和样条插值的类别，生成三次样条插值函数
def Cubic_Spline(xi, f, f_, kind):
    yi = f(xi)
    y_i = f_(xi)


    # 根据公式求出4个参数
    h = np.array([xi[k+1]- xi[k] for k in range(len(xi) - 1)],dtype = np.float32)
    u = np.array([h[k-1]/(h[k-1] + h[k]) for k in range(len(h) - 1)] + [1.0])
    lam = np.array([1.0] + [h[k]/(h[k-1] + h[k]) for k in range(len(h) - 1)])

    n = len(xi)
    d = np.zeros(n) + 1
    dyi = np.array([(yi[k+1] - yi[k])/h[k] for k in range(len(xi) - 1)])
    for j in range(1,n - 1):
        d[j] = 6 * (dyi[j] - dyi[j - 1])/(h[j-1] +h[j] )

    #根据参数判断是第几类样条插值 
    #第一类则使用题目给定的两端点一次导数值
    if kind == 'first':
        d[0] = 6/h[0] * (dyi[0] - 1)
        d[-1] = 6/h[-1] * (1 - dyi[-1])
     #第二类则使用原函数两端点一次导数值
    else:
        d[0] = 6/h[0] * (dyi[0] - y_i[0])
        d[-1] = 6/h[-1] * (y_i[-1] - dyi[-1])

    # 利用矩阵运算得出M
    A = np.mat(np.eye(n,n)) * 2.0
    for i in range(n-1):
        A[i + 1,i] = u[i]
        A[i,i+1] = lam[i]

    B = np.matrix(d)
    B = B.T

    M =  A.I * B
    M = np.array(M.T[0])[0]

    # 生成三次样条插值函数
    def Sx(x):
        result = 0.0
        # 利用k来定位x所处的区间
        k = 0
        while k < n and not (xi[k] <= x and xi[k+1] >= x):
            k += 1
        result = M[k] * (xi[k+1] - x)**3 / (6*h[k]) + M[k+1] * (x - xi[k])**3 / (6*h[k]) +\
            (yi[k] - (M[k]*h[k]**2)/6) * ((xi[k+1] - x) / h[k]) + \
            (yi[k+1] - (M[k+1]*h[k]**2)/6) * ((x - xi[k])/h[k])
        return result
    return Sx

第二类三次样条插值：

算法描述：

这里只是把两个端点的导数值换成原函数的导数值

本题条件：

f' (- 5) = 1 2

$f'(-5)=\frac12$

f' (5) = - 1 2

$f'(5)=-\frac12$

效果图如下：

这里写图片描述
可见在这种情况之下图像基本完全重合，拟合性很优越

总体平均误差值:

0.00402051060792

插值函数	误差值
Lagrange	0.285755407117
Hermite	0.350297537083
ChebyshevX	0.0486976071887
subLiner	0.0147779905191
subHermite	0.00132363164301
First Cubic Spline	0.0313803197288
Second Cubic Spline	0.00402051060792

从数值来看，各个积分公式的误差大小排序如下：

0.00132363164301 < <script type="math/tex" id="MathJax-Element-172"><</script> 0.00402051060792 < <script type="math/tex" id="MathJax-Element-173"><</script> 0.0147779905191 < <script type="math/tex" id="MathJax-Element-174"><</script> 0.0313803197288 < <script type="math/tex" id="MathJax-Element-175"><</script> 0.0486976071887 < <script type="math/tex" id="MathJax-Element-176"><</script> 0.285755407117 < <script type="math/tex" id="MathJax-Element-177"><</script> 0.350297537083
即：
subHermite < Second Cubic Spline < subLiner < First Cubic Spline < ChebyshevX < Lagrange < Hermite

根据表格中各个插值函数的平均误差值来看，效果最优的插值函数为分段埃尔米特插值函数

2.用以下数值积分公式计算积分

\int 1 - 1 1 1 + x 2 d x

$\int_{-1}^{1} \frac{1}{1+x^2}dx$

(1) 将区间 $[-1,1]$ 等分成 $20$ 份，用复合的梯形公式计算

算法描述：

把积分的区间 $[a, b]$ 分成 $N$ 份，分割出的每一个区间长度是一致的，然后就应用梯形公式：

\int_{a}^{b} f (x) d x \approx \frac{h}{2} \sum_{k = 0}^{n - 1} [f (x_{k}) + f (x_{k + 1})]

$\int_a^bf(x)\,dx\approx\frac h2\sum_{k=0}^{n-1}\left[f(x_k)+f(x_{k+1})\right]$

误差计算：

使用数学方法计算出

I = \int 1 - 1 1 1 + x 2 d x = π 2

$I=\int_{-1}^{1} \frac{1}{1+x^2}dx=\frac\pi2$
通过复合梯形公式计算出来的值为：

T = 1.56996299445

$T=1.56996299445$
误差：

| I - T | = 0.000833332341317

$\left|I-T\right|=0.000833332341317$

实现程序：

# 复合梯形公式
def Compound_trapezoid(x,h,f):
    result = 0.0
    n = len(x) - 1
    for k in range(n):
        result += f(x[k]) + f(x[k+1])
    result *= h/2
    return result

(2) 将区间 $[-1,1]$ 等分成 $10$ 份，用复合的辛普森公式计算

算法描述：

将区间 $[a,b]$ 分成 $n$ 等份，在每个子区间 $[x_k,x_{k+1}]$ 上采用辛普森公式：

\int b a f (x) d x \approx = h 6 \sum k = 0 n - 1 [f (x k) + 4 f (x k + 1 / 2) + f (x k + 1)]

$\int_a^bf(x)\,dx\approx=\frac h6\sum_{k=0}^{n-1}\left[f(x_k)+4f\left(x_{k+1/2}\right)+f(x_{k+1})\right]$

误差计算：

通过复合的辛普森公式计算出来的结果:

S = 1.57079630697

$S = 1.57079630697$

误差：

| I - S | = 0.0000000198252885

$\left|I-S\right|=0.0000000198252885$

实现程序：

# 复合辛普森公式
def Compound_Simpson(x,h,f):
    result = 0.0
    n = len(x) - 1
    for k in range(n):
        result += f(x[k]) + 4*f(x[k] + h/2) + f(x[k+1])
    result *= h/6
    return result

(3) 将区间 $[-1,1]$ 等分成 $10$ 份，用复合的两点高斯公式计算

算法描述：

在高斯求积公式中，若取权函数 $\rho(x)=1$ ，区间为 $[-1,1]$ ，则得到公式：

\int 1 - 1 f (x) d x \approx \sum k = 0 n A k f (x k)

$\int_{-1}^1f(x)\,dx\approx\sum_{k=0}^nA_kf(x_k)$
当积分区间不是

[−1,1] [ − 1 , 1 ] $[-1,1]$ ，而是一般的积分区间

[a,b] [ a , b ] $[a,b]$ 时，需做以下变换：

x = b - a 2 t + a + b 2

$x=\frac{b-a}2t+\frac{a+b}2$
可将

[a,b] [ a , b ] $[a,b]$ 化为

[−1,1] [ − 1 , 1 ] $[-1,1]$ ，这时

\int b a f (x) d x = b - a 2 \int 1 - 1 f (b - a 2 t + a + b 2) d t

$\int_a^bf(x)\operatorname dx=\frac{b-a}2\int_{-1}^1f\left(\frac{b-a}2t+\frac{a+b}2\right)\operatorname dt$
将

[−1,1] [ − 1 , 1 ] $[-1,1]$ 划分成

10 10 $10$ 个等份之后，将每一部分通过区间变化成

[−1,1] [ − 1 , 1 ] $[-1,1]$ ，接着将每一部分积分出来的值相加即可得到积分值

误差计算：

通过公式计算得:

G 2 = 1.5707964042

$G_2=1.5707964042$
相对误差

| I - G 2 | = 0.0000000774061313

$\left|I-G_2\right|=0.0000000774061313$

实现程序：

#查表得出各个情况的Xk和Ak
Xk = [
    [0.0000000],
    [-0.5773503, 0.5773503],
    [-0.7745967, 0.0000000, 0.7745967],
    [-0.8611363, -0.3399810, 0.3399810, 0.8611363],
    [-0.9061798, -0.5384693, 0.0000000, 0.5384693, 0.9061798],
    [-0.9324695, -0.6612094, -0.2386192, 0.2386192, 0.6612094, 0.9324695]
]

Ak = [
    [2.0000000],
    [1.0000000, 1.0000000],
    [0.5555556, 0.8888889, 0.5555556],
    [0.3478548, 0.6521452, 0.6521452, 0.3478548],
    [0.2369269, 0.4786287, 0.5688889, 0.4786287, 0.2369269],
    [0.1713245, 0.3607616, 0.4679139, 0.4679139, 0.3607616, 0.1713245]
]

# 复合高斯公式
def Compound_Goss(point, f, begin=-1.0, end=1.0, n=1):
    # 根据传入参数确定是哪种高斯公式，选择对应的Xk和Ak值
    x = Xk[point]
    A = Ak[point]

    # 计算间隔h
    h = (end - begin) / n
    # print h
    I = np.zeros(n,dtype=np.float64)
    for i in range(n):
        a = begin + i*h
        b = begin + (i+1)*h
        # print a,b
        I[i] = sum([ f(change(a,b,x[k])) * A[k] for k in range(len(x)) ])
        I[i] *= (b-a)/2
    return sum(I)

(4) 将区间 $[-1,1]$ 等分成 $6$ 份，用复合的三点高斯公式计算

算法描述：

与复合的两点高斯公式计算方法类似
将 $[-1,1]$ 划分成 $6$ 个等份之后，将每一部分通过区间变化成 $[-1,1]$ ，接着将每一部分积分出来的值相加即可得到积分值

误差计算：

通过公式计算得:

G 3 = 1.57079632759

$G_3=1.57079632759$
相对误差

| I - G 3 | = 0.0000000007934180

$\left|I-G_3\right|=0.0000000007934180$

(5) 将区间 $[-1,1]$ 等分成 $4$ 份，用复合的五点高斯公式计算

算法描述：

与复合的两点高斯公式计算方法类似
将 $[-1,1]$ 划分成 $4$ 个等份之后，将每一部分通过区间变化成 $[-1,1]$ ，接着将每一部分积分出来的值相加即可得到积分值

误差计算：

通过公式计算得:

G 5 = 1.57079632614

$G_5=1.57079632614$
相对误差

| I - G 5 | = 0.0000000006552323

$\left|I-G_5\right|=0.0000000006552323$

误差比较与评价

函数	积分值	误差
复合的梯形公式	1.56996299445	0.0008333323413172
复合的辛普森公式	1.57079630697	0.0000000198252885
复合的两点高斯公式	1.5707964042	0.0000000774061313
复合的三点高斯公式	1.57079632759	0.0000000007934180
复合的五点高斯公式	1.57079632614	0.0000000006552323

从数值来看，各个积分公式的误差大小排序如下：

0.0000000006552323 < <script type="math/tex" id="MathJax-Element-225"><</script> 0.0000000007934180 < <script type="math/tex" id="MathJax-Element-226"><</script> 0.0000000198252885 < <script type="math/tex" id="MathJax-Element-227"><</script> 0.0000000774061313 < <script type="math/tex" id="MathJax-Element-228"><</script> 0.0008333323413172
即：
复合的五点高斯公式 < <script type="math/tex" id="MathJax-Element-229"><</script> 复合的三点高斯公式 < <script type="math/tex" id="MathJax-Element-230"><</script> 复合的辛普森公式 < <script type="math/tex" id="MathJax-Element-231"><</script> 复合的两点高斯公式 < <script type="math/tex" id="MathJax-Element-232"><</script> 复合的梯形公式

可见复合的五点高斯公式误差值最小，同时复合的五点高斯公式和复合的三点高斯公式都已经把误差控制在 $10^{-10}$ 数量级，误差已经非常小了。
而复合的梯形公式的误差在 $10^{-4}$ ，误差相对来说就非常大

Angra-Mainyu

关注

26
点赞
踩
207

收藏

觉得还不错? 一键收藏
3
评论
数值分析实习作业（各种插值函数与积分公式的python代码实现）

数值分析实习作业院系学号姓名数据科学与计算机学院 16339053 张志宏1.对函数f(x)f(x)f(x)进行插值f(x)=11+x2f(x)=11+x2f(x)=\frac{1}{1+x^{2}}(1)令插值节点为等距节点−5,−4,−3,−2,−1,0,1,2,3,4,5−5,−4,−3,−2,−1,0,1,2,3,4,...
复制链接

扫一扫