数值分析常见习题解答

1．已知下列数值表，求符合表值的插值多项式，并给出插值余项的表达式。

xix_ixi	000	111	222
yiy_iyi	222	111	222
yi′y_i^{'}yi′	−2-2−2	−1-1−1
yi′′y_i^{''}yi′′	−10-10−10

x_i

y_i

y_i^{'}

-2

-1

y_i^{''}

-10

解：采用牛顿插值：

P_2(x)=f(x_0)+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)

=x^2-2x+2

由题目条件可得符合该表的插值多项式可设为：

P_5(x)=P_2(x)+(ax^2+bx+c)x(x-1)(x-2)

代入以下条件，有：

KaTeX parse error: Unknown column alignment: 1 at position 28: … begin{array}{1̲} P…

我们可以得到：

P_5(x)=4x^5-15x^4+17x^3-5x^2-2x+2

R(f)=frac{f^{(6)}(xi)}{6!}x^3(x-1)^2(x-2),式中xi位于x_0,x_1x_2和x所界定的范围内

解：

由给定条件，可确定不超过三次的插值多项式，其形式为：

P(x)=f(x_0)+f[x_0,x_1](x-x_0)+A(x-x_0)(x-x_1),代入P^{'}(x_1)=f^{'}(x_1),得

A=frac{f^{'}(x_0)-f[x_0,x_1]}{x_0-x_1}

代入P_(x)，得:

P_(x)=-frac{(x-x_1)(x-2x_0+x_1)}{(x_1-x_0)^2}f(x_0)+frac{(x-x_0)(x-x_1)}{x_0-x_1}f^{'}(x_0)+frac{(x-x_0)^2}{(x_1-x_0)^2}f(x_1)

为了求出余项R(x)=f(x)-P(x),设R(x)=f(x)-P(x)=k(x)(x-x_0)^2(x-x_1),其中k(x)为待定常数

构造 varphi(t)=f(t)-P(t)-k(x)(t-x_0)^2(t-x_1),thereforevarphi(t)在(x_0,x_1)内有四个零点，根据罗尔定理结论，可得\varphi^{(3)}(t)在(x_0,x_1)内至少有一个零点xi,使得：

varphi^{(3)}(xi)=f^{(3)}(xi)-3!k(x)=0

therefore k(x)=frac{f^{(3)}(xi)}{6},R(x)=frac{1}{6}(x-x_0)^2(x-x_1)f^{(3)}(xi),x_0<xi<x_1,证毕

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-z4NG02Jk-1587562372068)(https://s1.ax1x.com/2020/04/21/J8zKTU.png)]

解：

法一：构造正交多项式varphi_0(x),varphi_1(x),varphi_2(x)求解

qquad varphi_0(x)=1\ qquad a_0=frac{(f,varphi_0)}{(varphi_0,varphi_0)}=14.5\ qquad a_1=frac{(xvarphi_0,varphi_0)}{(varphi_0,varphi_0)}=2.5\ qquadvarphi_1(x)=(x-a_1)varphi_0(x)\qquad=x-2.5 \ quadbeta_1=frac{(varphi_1(x),varphi_1(x))}{(varphi_1(x),varphi_1(x))}\ qquadvarphi_2(x)=(x-a_1)varphi_1(x)-beta_1varphi_0(x)\=x^2-5x+5\ a_2=frac{(xvarphi_1,varphi_1)}{(varphi_1,varphi_1)}=0.5\

therefore y=a_0varphi_0(x)+a_1varphi_1(x)+a_2varphi_2(x)=frac{1}{2}x^2+frac{49}{10}x-frac{3}{2}

法二：取phi=span(1,x,x^2),求解线性方程组，也可得到相同答案，不再赘述

4.求f(x)=x^4+3x^3-1在区间[0,1]上的三次最佳一致逼近多项式

解:

做变换;x=frac{t+1}{2},有:f(t)=frac{1}{16}(1+t)^4+frac{3}{8}(t+1)^3-1

由首项为1的切比雪夫多项式性质可得：

f(x)-P_3(x)=frac{1}{16}cdotfrac{1}{8}T_4(t)

therefore P_3(x)=f(x)-frac{1}{128}T_4(t)=5x^3-frac{5}{4}x^2+frac{1}{4}x-frac{129}{128},xin[0,1]

解:

由代数精度定义可知，分别令f(x)=1,x-x_0,(x-x_0)^2,(x-x_0)^3,得：

KaTeX parse error: Unknown column alignment: 1 at position 28: … begin{array}{1̲} A…

解得：

KaTeX parse error: Unknown column alignment: 1 at position 28: … begin{array}{1̲} A…

同时有：

int_{x_0}^{x_1}(x-x_0)(x-x_0)^4dxneq h ^2(0+Bh^4)+h^3(0+4h^3D)

所以该数值积分具有三次代数精度，余项为:

R(f)=int_{x_0}^{x_1}frac{f^{(4)}(xi)}{4!}(x-x_0)^3(x-x_1)^2=f^{(4)}(eta)frac{h^2}{4!cdot60}=frac{h^2}{1440}f^{(4)}(eta),xin[x_0,x_1],xiin(x_0,x_1),etain(x_0,x_1)

证明：

对于插值基函数，有l_i(x_j)=delta_{ij}

thereforeint_a^b rho(x)l_k^2(x)dx=sum_{i=0}^nA_il_k^2(x_i)=A_k

sum_{k=0}^nint_a^b rho(x)l_k^2(x)dx=sum_{k=0}^nA_k=int_a^brho(x)dx

7．用龙贝格方法计算 int_{0.3}^{0.8}x^3+frac{sinx}{x}, 使其误差不超过 0.1*10^{-4}:

解：编写Python程序：

import math

'''
给定一个函数，如：f(x)= x^(3/2），和积分上下限a,b，用机械求积Romberg公式求积分。

'''
import numpy as np


def func(x):
    return x**2+math.sin(x)/x

class Romberg:
    def __init__(self, integ_dowlimit, integ_uplimit):
        '''
        初始化积分上限integ_uplimit和积分下限integ_dowlimit
        输入一个函数，输出函数在积分上下限的积分

        '''
        self.integ_uplimit = integ_uplimit
        self.integ_dowlimit = integ_dowlimit



    def calc(self):
        '''
        计算Richardson外推算法的四个序列

        '''
        t_seq1 = np.zeros(5, 'f')
        s_seq2 = np.zeros(4, 'f')
        c_seq3 = np.zeros(3, 'f')
        r_seq4 = np.zeros(2, 'f')
        # 循环生成hm间距序列
        hm = [(self.integ_uplimit - self.integ_dowlimit) / (2 ** i) for i in range(0,5)]
        print(hm)
        # 循环生成t_seq1
        fa = func(self.integ_dowlimit)
        fb = func(self.integ_uplimit)

        t0 = (1 / 2) * (self.integ_uplimit - self.integ_dowlimit) * (fa+fb)
        t_seq1[0] = t0

        for i in range(1, 5):
            sum = 0
            # 多出来的点的累加和
            for each in range(1, 2**i,2):
                sum =sum + hm[i]*func( self.integ_dowlimit+each * hm[i])#计算两项值
            temp1 = 1 / 2 * t_seq1[i - 1]
            temp2 =sum
            temp =  temp1 + temp2
            # 求t_seql的1-4位
            t_seq1[i] = temp
        print('T序列：'+ str(list(t_seq1)))
        # 循环生成s_seq2
        s_seq2 = [round((4 * t_seq1[i + 1] - t_seq1[i]) / 3,6) for i in range(0, 4)]
        print('S序列：' + str(list(s_seq2)))
        # 循环生成c_seq3
        c_seq3 = [round((4 ** 2 * s_seq2[i + 1] - s_seq2[i]) / (4 ** 2 - 1),6) for i in range(0, 3)]
        print('C序列：' + str(list(c_seq3)))
        # 循环生成r_seq4
        r_seq4 = [round((4 ** 3 * c_seq3[i + 1] - c_seq3[i]) / (4 ** 3 - 1),6) for i in range(0, 2)]
        print('R序列：' + str(list(r_seq4)))
        r_seq5 = [round((4 ** 4 * r_seq4[i + 1] - r_seq4[i]) / (4 ** 4 - 1),6) for i in range(0, 1)]
        print('A序列：' + str(list(r_seq5)))
        return 'end'


rom = Romberg(0.3, 0.8)
print(rom.calc())

运行结果：

、

最终我们得到，int_{0.3}^{0.8}frac{x^3+sinx}{x}=0.635258

1),6) for i in range(0, 1)] print(‘A序列：’ + str(list(r_seq5))) return ‘end’

rom = Romberg(0.3, 0.8) print(rom.calc())

$运行结果：$

[外链图片转存中...(img-l6ySF4Lu-1587562372076)]、

$最终我们得到，int_{0.3}^{0.8}frac{x^3+sinx}{x}=0.635258$