-数值分析-

一、引论

（一）数值分析的对象、作用及特点

1.1.1 数值分析的对象

解决问题的步骤通常如下：

1、建立数学模型
2、由模型给出数学计算方法
3、编制程序计算结果

第一点是应用数学的任务，第二、三点是计算数学的任务，也是数值分析的研究对象。

1.1.2 数值分析的特点

数值问题：输入数据与输出数据之间函数关系的一个确定而无歧义的描述，输入与输出数据用有限维数据表示

计算的基本单位是算法元

算法元=算子+输入元+输出元

数值分析是研究数值问题的算法，它的特点：

• 面向计算机
• 有可靠的理论分析
• 有较好的计算复杂性（time space）
• 要有数值实验

（二）数值计算的误差

1.2.1 误差来源及分类

• 模型误差：数学模型与实际问题之间出现的误差
• 观测误差：使用仪器观测数据时产生的误差
• 截断误差/方法误差：利用数值方法逼近真实解时的误差
• 舍入误差：由计算机字长存储所导致的误差

1.2.2 误差及有效数字

设x为准确值，x为近似值
绝对误差：e=x-x*
绝对误差限：|x-x*|≤ε*
相对误差：e*/x，但由于x无法知道，用e_r^*=e*/x*
相对误差限：ε_r^*=ε*/x*
注意：绝对误差和绝对误差限是有量纲的，相对误差和相对误差限没有量纲

一般地，凡是由四舍五入方法得到的近似值，其绝对误差限等于近似值末位的半个单位

有效数字：若近似值x的误差限为某一位的半个单位，该位到x的第一位非零数字有n位，则x*有n位有效数字，表示为：

x*=±10^m*(a₁+a₂*10^-2+…+a_n*10^-n+1)，其中a_i为0-9的数字，a₁≠0，且|x*-x|≤1/2*10^m-n+1

有效数字与相对误差限的关系：
对于x*，若其有n位有效数字，则相对误差限ε_r^*≤(1/2a₁)*10^-n+1
若相对误差限ε_r^*≤(1/2(a₁+1))*10^-n+1，则x*的有效位数为n位

1.2.3 数值运算的误差估计（略）

设精确值x₁和x₂的近似值为x₁^*和x₂^*，其+、-、*、/满足以下不等式：
ε(x₁^*±x₂^*)≤ε(x₁^*)+ε(x₂^*);
ε(x₁^*x₂^*)≤|x₁^*|ε(x₂^*)+ε(x₁^*)|x₂^*|;
ε(x₁^*/x₂^*)≤(|x₁^*|ε(x₂^*)+ε(x₁^*)|x₂^*|)/|x₂^*|²;
对于单元函数：
ε(f(x^*))≈|f’(x*)|ε(x^*)
对于多元函数：
ε(A^*)≈[Σ_k=1ⁿ|(∂f/∂x_k)*|]ε(x^*)

（三）误差定性分析与避免误差危害

1.3.1 算法数值稳定性

算法的数值稳定性：一个算法的输入数据如果有误差，而在计算舍入误差不增长，则说该算法是稳定的，否则反之

例：计算I_n=∫₀¹ xⁿe^x-1dx并估计误差
使用分部积分可以得到：

• I_n=1-nI_n-1
• I₀=∫₀¹e^x-1dx=1-e^-1

当我们用4为有效数字表示I₀时会发现I₈=-0.7280误差很大，并且与结果正负性相悖，这是因为，通过递推可以得到e_n=(-1)ⁿn!e₀，递推n次以后的误差是一开始的阶乘次级，因此该算法是不稳定的

1.3.2 病态问题与条件数

病态问题：对于一个数值问题本身，若其输入数据有微小误差却导致输出数据的巨大误差，则说这是病态问题。

条件数：函数值f(x*)的相对误差和x*的相对误差即为条件数

|(f(x)-f(x*))/f(x)|/|Δx/x|≈|xf’(x)/f(x)|=C_p

注意：病态问题是数值问题本身固有的，其条件数越大(10≤C_p即为病态问题)，病态问题越严重

1.3.3 避免误差危害

1. 避免相近两数相减：相近两数相减会丢失有效位数
   常用方法：

(1)替代法
logx₁-logx₂=logx₁/x₂
1-cosx=2sin²x/2
√(x+1)-√(x)=1/(√(x+1)+√(x))
(2)采用双倍计算机字长计算

2. 防止大数吃掉小数：同量级数字会因为字长限制失去小数的精度从而导致失真，例如A=23151671+0.2在八位计算机对阶会丢失0.2
3. 避免用绝对值较小的数作除数：溢出且绝对误差很大
4. 注意简化计算程序，减少计算次数
5. 选用数值稳定性好的算法

二、插值法

（一）引言（略）

2.1.2 多项式插值

现实问题中y=f(x)表示某种内在规律的关系，但f(x)无法准确得知，仅知道y=f(x)在[a,b]上有定义，且已知a≤x₀＜x₁＜···＜x_n≤b的值为y₀,y₁,···,y_n，利用一个简单函数P(x)满足P(x_i)=y_i成立，这种求取P(x)的方法叫做插值法。
其中，P(x)为f(x)的插值函数，[a,b]为插值区间，x₀,x₁,···,x_n为插值节点。根据P(x)的类型可以分为多项式插值，分段插值，三角插值等
从几何上看，插值法就是求取一个函数拟合真实曲线的方法

满足P(x_i)=y_i的插值多项式存在且唯一：

存在性：
对于多项式插值可以得到以下关于a₀,a₁,···,a_n的n+1元方程组：

其系数矩阵：
由于x_i之间互不相同，所以D≠0，系数矩阵不为0即可证明P(x)一定有解

唯一性：（反证法）
假设存在两个插值函数p₁(x),p₂(x)满足条件，那么对于每个x_i都有p₁(x_i)-P(x_i)=0，此时p₁(x)存在n+1个零点，这和n次多项式只能有n个零点相悖，因此p₁(x)=p₂(x)

（二）拉格朗日插值

2.2.1 线性插值和抛物线插值

线性插值：
已知两点(x₀,y₀),(x₁,y₁)拟合曲线f(x):

两点式：L(x)=y₀*(x-x₁)/(x₀-x₁)+y₁*(x-x₀)/(x₁-x₀)

我们把(x-x₁)/(x₀-x₁)用l₀(x)表示，(x-x₀)/(x₁-x₀)用l₁(x)表示，即

L(x)=y₀l₀(x)+y₁l₁(x)
把l_i(x)称为线性插值基函数，其满足：

l₀(x₀)=1,l₀(x₁)=0
l₁(x₀)=0,l₁(x₁)=1

抛物线插值：
已知三点(x₀,y₀),(x₁,y₁),(x₂,y₂)拟合曲线:

利用线性插值中的基函数方法：
L(x)=y₀l₀(x)+y₁l₁(x)+y₂l₂(x)
其中l_i(x)满足：

• l₀(x₀)=1,l₀(x₁)=0,l₀(x₂)=0
• l₁(x₀)=0,l₁(x₁)=1,l₁(x₂)=0
• l₂(x₀)=0,l₂(x₁)=0,l₂(x₂)=1

由上知l₀(x)有两个零点，因此可以设l₀(x)=A(x-x₁)(x-x₂)，再由l₀(x₀)可知A=1/(x₀-x₁)(x₀-x₂)，得l₀(x)=(x-x₁)(x-x₂)/(x₀-x₁)(x₀-x₂)，同理 l₁(x)，l₂(x)

2.2.2 拉格朗日多项式插值

把线性插值和抛物线插值推广到n次多项式的情况即可得到n次多项式插值：

L_n(x)=Σ_k=0ⁿy_kl_k(x)
l_k(x)=(Σ_i=0^k-1(x-x_i)*Σ_i=k+1ⁿ(x-x_i))/((Σ_i=0^k-1(x_k-x_i)**Σ_i=k+1ⁿ(x_k-x_i)))
即分子为x-任何点，除了x_k的点，分母为x_k-任何点，除了x_k的点

若记ω_n+1(x)=(x-x₀)(x-x₁)···(x-x_n)
则ω’_n+1(x_k)=(x_k-x₀)···(x_k-x_k-1)(x_k-x_k+1)···(x_k-x_n)
L_n(x)=Σ_k=0ⁿω_n+1(x)/((x-x_k)ω’_n+1(x_k))

2.2.3 插值余项与误差估计

在使用L_n(x)近似f(x)时会产生误差，若f⁽ⁿ⁾(x)在[a,b]连续且f⁽ⁿ⁺¹⁾(x)在(a,b)存在（前提），那么其截断误差R_n(x)=f(x)-L_n(x)称为插值多项式的余项，也是误差，且

R_n(x)=f(x)-L_n(x)=(f⁽ⁿ⁺¹⁾(ξ)/(n+1)!)ω_n+1(x)，ξ∈(a,b)

证明如下：

在x取x₀,x₁,···,x_n时，L_n(x)是准确的，即R_n(x)=0，除此之外的x的误差另外计算，因此设
R_n(x)=K(x)(x-x₀)(x-x₁)···(x-x_n)=K(x)ω_n+1(x),K(x)是与x有关的函数
现在要求出K(x)，确定x，新作一个函数：
φ(t)=f(t)-L_n(t)-K(x)(x-x₀)(x-x₁)···(x-x_n)
由前提可知φ⁽ⁿ⁾(t)在[a,b]连续，并且φ⁽ⁿ⁺¹⁾(t)在[a,b]存在，可以知道φ(t)在x₀,x₁,···,x_n和x时等于0，即存在n+2个零点，根据罗尔定理（见下）可知，φ‘(t)在φ(t)两个零点间至少存在一个零点，因此φ’(t)至少存在n+1个零点，同理φ‘’(t)至少存在n个零点，依此类推，φ⁽ⁿ⁺¹⁾(t)至少存在一个零点，则：
令φ⁽ⁿ⁺¹⁾(ξ)=f⁽ⁿ⁺¹⁾(ξ)-(n+1)!K(x)=0
K(x)=f⁽ⁿ⁺¹⁾(ξ)/(n+1)!
得证

罗尔定理：如果 R 上的函数 f (x) 满足以下条件：（1）在 闭区间 [a,b] 上 连续 ，（2）在 开区间 (a,b) 内 可导 ，（3）f (a)=f (b)，则至少存在一个 ξ∈ (a,b)，使得 f’ (ξ)=0。

由于ξ值是无法知道的，因此我们可以通过放缩来逼近截断误差限，令max|f⁽ⁿ⁺¹⁾(x)|=M_n+1，x∈[a,b]

|R_n(x)|≤f(x)-L_n(x)=(M/(n+1)!)|ω_n+1(x)|

当n=1时，R₁(x)=1/2 * f’‘(ξ)(x-x₀)(x-x₁)
当n=2时，R₂(x)=1/6 * f’‘'(ξ)(x-x₀)(x-x₁)(x-x₂)
当f(x)=x^k,k≤n时，f⁽ⁿ⁺¹⁾=0，R_n(x)=x^k-Σ_i=0ⁿx_i^kl_i(x)=0，x^k=Σ_i=0ⁿx_i^kl_i(x)

（三）均差与牛顿插值多项式

拉格朗日插值法的局限性：当节点的个数确定且不变时，拉格朗日多项式是一个很好的方法，但是一旦数据节点发生变化或节点个数增加，插值多项式就需要重新计算。我们需要一个新的方法使插值多项式逐次生成。

2.3.2 均差及其性质

线性插值除了用两点式，还可以用点斜式来表示：

P₁(x)=y₀+y₁*(y₁-y₀)/(x₁-x₀)

我们把(y₁-y₀)/(x₁-x₀)记作a₁，推广到抛物线插值：

P₂(x)=y₀+y₁*a₁+y₂a₂
它满足：
P₂(x₀)=y₀，P₂(x₁)=y₁，P₂(x₂)=y₂
解得：
a₁=(y₁-y₀)/(x₁-x₀)
a₂=(((y₂-y₀)/(x₂-x₀))-((y₁-y₀)/(x₁-x₀)))/(x₂-x₁)

为了便于定义，给出均差（茶商）的定义：

称f[x₀,x_k]=(f(x_k)-f(x₀))/(x_k-x₀)为f(x)关于x_k,x₀的一阶差商
称f[x₀,x₁,x_k]=(f[x₀,x_k]-f[x₀,x₁])/(x_k-x₁)为f(x)关于x_k,x₁,x₀的二阶差商
称f[x₀,x₁,···,x_k]=(f[x₀,···,x_k-2,x_k]-f[x₀,···,x_k-2,x_k-1])/(x_k-x_k-1)为f(x)关于x_k,x₁,x₀的k阶差商

均差的三个性质

1. k阶均差可以表示成f(x₀),f(x₁),···,f(x_k)的线性组合
   f[x₀,x₁,···,x_k]=Σ_j=0^k(f(x_j)/(x_j-x₀)···(x_j-x_j-1)(x_j-x_j+1)···(x_j-x_k))
   该性质也表明均差与节点次序无关（均差的对称性）
2. f[x₀,x₁,···,x_k]=(f[x₁,x₂,···,x_k]-f[x₀,x₁,···,x_k-1])/(x_k-x₀)
3. 若f(x)存在n阶导数，则f[x₀,x₁,···,x_k]=f^(n)(ξ)/n!

2.3.3 牛顿插值多项式

我们把P₀(x)设为y₀，那么
P₀(x)=y₀
P₁(x)=P₀(x)+fx₀,x₁
P₂(x)=P₁(x)+fx₀,x₁,x₂(x-x₁)
···
P_n(x)=P_n-1(x)+fx₀,x₁,···,x_k(x-x₁)···(x-x_n-1)
其余项和拉格朗日插值余项是一致的
R_n(x)=f(x)-P_n(x)=f[x₀,x₁,···,x_k]ω_n+1(x)

牛顿插值多项式相较于拉格朗日多项式更具有一般性，它允许导数不存在及离散点的情况

（五）分段低次插值

2.5.1 高次插值的病态性质

龙格（Runge）现象：当使用高次函数拟合实际曲线时，当n趋于无穷时，高次函数P_n(x)不一定收敛于L(x)（或截断误差很大）

对于L_n(x)=1/(1+x²)使用插值方法计算拟合函数，可以看到图像与实际图像的拟合程度很差，这便是龙格现象（可以理解为过拟合数据点或模型复杂度太高）其图像如图：

2.5.2 分段线性插值

为解决龙格现象，我们把数据进行分段，利用低次的折线段来逼近实际函数曲线，即

对于在[a,b]区间上的x₀,x₁,···,x_n，对每个[x_k,x_k+1]，用
I_h(x)=y_k(x-x_k+1)/(x_k-x_k+1)+y_k+1(x-x_k)/(x_k+1-x_k)，k=0,1,···,n-1
来代替，得到的 I_h(x)称为分段线性插值函数

其误差估计用插值余项可得：

设h_k=x_k+1-x_k，h=max|h_k|
max|f(x)-I_h(x)|≤(M₂/2)*max|(x-x_k)(x-x_k-1)|，
放缩后得max|f(x)-I_h(x)|≤M₂h²/8
这里x_k≤x≤x_k+1，M₂=max|f’'(x)|

（六）三次样条插值

2.6.1 三次样条函数

分段线性插值具有收敛性，但其光滑性差，对于某些模型如船体模型，机翼样形都需要较高的光滑度，要求有二阶连续导数，因此提出三次样条函数来拟合此类曲线

三次样条函数：对于[a,b]区间上的x₀,x₁,···,x_n，其每一区间[x_k,x_k+1]是一个三次函数，所构成的属于C²[a,b]的S(x)是三次样条函数，若对于每个y_i=f(x_i)都有y_i=S(x_i)，则S(x)为三次样条插值函数。

对于这样一个三次样条插值函数S(x)，其在[a,b]上有n个区间，每一个区间用一个三次多项式ax³+bx²+cx+d来表示，即含义四个未知数，在整个区间上共有4n个未知数，一个未知数需要一个方程（一个条件），因此需要4n个条件，在S(x)中，其二阶导数光滑，且满足插值，则有：
S(x_j-0)=S(x_j+0)
S’(x_j-0)=S’(x_j+0)
S’‘(x_j-0)=S’'(x_j+0)
共有n-1个非边界点，于是有3n-3个条件，又
S(x_i)=y_i
有n+1个条件，此时有4n-2个条件，还欠缺两个条件，往往会以下列三种形式给出：

（1）已知两端导数值：f’(x₀)=f₀’,f’(x_n)=f_n’
（2）已知两端二阶导数值：f‘‘(x₀)=f₀’’,f’‘(x_n)=f_n’’
（3）函数呈以x_n-x₀为周期的形式：

S(x_j-0)=S(x_j+0)
S’(x_j-0)=S’(x_j+0)
S’‘(x_j-0)=S’'(x_j+0)

2.6.2 样条插值函数的建立

三转角方法：
（当给出f‘(x₀)和f’(x_n)时更方便）设S’(x_j)=m_j，使用分段三次埃尔米特插值，由插值条件可得

S(x)=Σ_j=0ⁿ[y_jα_j(x)+m_jβ_j(x)]

α_j(x)和β_j(x)被称为插值基函数，由于各项除了m都是已知的，求解S(x)的问题就转变为求解m_j的问题

在分段三次埃尔米特插值中，把区间分成n-1个区间，从[x₀,x₁],[x₁,x₂]到[x_n-1,x_n]，>对x∈[x_i-1,x_i]得到
S(x)=[(x-x_i)²(x_i+2x-3x_i-1)y_i-1+(x-x_i-1)²(3x_i-2x-3x_i)y_i]/h_i³+[(x-x_i)²(x-x_i-1)m_i-1+(x-x_i-1)²(x-x_i)m_i]/h_i²，其中h_i=x_i-x_i-1
确定m_i，i=1,2,···,n-1，即可确定S(x)
利用S’‘(x_i-0)=S’'(x_i+0)即可列出关于m_i的方程
对S(x)求导两次得：
S‘’(x)=6(2x-x_i-1-x_i)(y_i-1-y_i)/h_i³+2[(3x-x_i-1-2x_i)m_i-1+(3x-2x_i-1-x_i)m_i]/h_i²
由前一项和后一项可得
λ_im_i-1+2m_i+μ_im_i+1=g_i，i=1,2,…,n-1
其中，λ_i=h_i+1/(h_i+h_i+1),μ_i=1-λ_i,g_i=3(λ_if[x_i,x_i]+μ_if[x_i,x_i+1])
当题目给出m₀和m_n时方程可求解

三弯矩方法：
（当给出S’‘(x₀)和S’‘(x_n)时更方便） 设S’'(x_i)由于S(x)是三次可微函数，其二阶导数为线性形式，可以设为

S’‘(x)=M_i-1(x-x_i)/(x_i-1-x_i)+M_i(x-x_i-`)/(x_i-x_i-1)
此时对S‘’(x)积分两次，就会得到两个常数项，再根据S(x_i-1)=y_i-1和S(x_i)=y_i确定积分常数，得：
S(x)=[(x_i-x)³M_i-1+(x-x_i-1)³M_i]/6h_i+(y_i-1/h_i-h_iM_i-1/6)(x_i-x)+(y_i/h_i-h_iM_i/6)(x-x_i-1)，其中h_i=x_i-x_i-1
同理，确定M_i即可确定S(x)，利用S’‘(x_i-0)=S’'(x_i+0)即可列出关于M_i的方程，可得
μ_iM_i-1+2M_i+λ_iM_i+1=d_i
λ和μ同上，d_i=6f[x_i-1,x_i,x_i+1]

三、函数逼近和快速傅里叶变换

（一）函数逼近的基本概念

在前面的插值问题中，我们通过精确点来确定一个插值多项式来模拟原来的函数，但实际中常常点是不精确的，我们通过用简单函数逼近已知函数的方法来使得f(x)和已知函数φ(x)的误差达到最小，这就是函数逼近问题

3.1.1 函数逼近与函数空间

在各种集合中引入的某些不同的确定关系称之为空间：

n维向量空间：由所有n维向量构成的集合按向量加法和数乘向量的乘法构成实数域上的线性空间
多项式空间：由次数不超过n的实系数多项式全体按多项式加法与数乘多项式规则构成的数域R上的一个线性空间，用H_n表示

线性代数相关知识：
对于x₁,x₂,···,x_n，如果存在不全为0的数α₁,α₂,···,α_n使得α₁x₁+α₂x₂+···+α_nx_n=0，则称x₁,x₂,···,x_n线性相关，否则线性无关，若任意x可以用x₁,x₂,···,x_n表示，x=α₁x₁+α₂x₂+···+α_nx_n，则x₁,x₂,···,x_n称为一组基，α₁,α₂,···,α_n称为在该基下的坐标

魏尔斯特拉斯定理：

设f(x)∈C[a,b]，则对于任何ε＞0，总存在一个代数多项式p(x)，使得max_a≤x≤b|f(x)-p(x)|＜ε在[a,b]上一致成立

3.1.2 范数与赋范线性空间

在上述定理中，|f(x)-p(x)|是简单函数与实际函数的距离，为了对线性空间中不同元素大小进行度量，使用范数来表示：

设S为线性空间，x∈S，若存在唯一实数|| ||，满足
（1）正当性：||x||≥0，当且仅当x=0时||x||=0
（2）齐次性：||αx||=|α| ||x||
（3）三角不等式：||x+y||≤||x||+||y||
则称|| ||为S上的范数，S与|| ||一起称为赋范线性空间

向量范数：

对于在Rⁿ上的向量x=(x₁,x₂,···,x_n)^T，
1-范数：||x||₁=Σ_i=1ⁿ|x_i|
2-范数：||x||₂=(Σ_i=1ⁿx_i²)^1/2
∞-范数：||x||_∞=max_1≤i≤n|x_i|

函数范数：

对于在连续空间C[a,b]上的，f∈C[a,b]，
1-范数：||f||₁=∫_a^b|f(x)|dx
2-范数：||f||₂=(∫_a^bf²(x)dx)^1/2
∞-范数：||f||_∞=max_a≤x≤b|f(x)|

3.1.3 内积与内积空间

线性代数中的内积：
(x,y)=x₁y₁+x₂y₂+···+x_ny_n
将其推广到一般的线性空间

设X为数域K上的线性空间，对于任意u,v∈X，有K中一个数与之对应，且满足：
（1）(u,v)=$\overline{\text{(v,u)}}$
（2）(αu,v)=α(u,v)
（3）(u+v,w)=(u,w)+(v,w)
（4）(u,u)≥0,当且仅当u=0时(u,u)=0
称(u,v)为X上u与v的内积，$\overline{\text{(v,u)}}$为(u,v)的共轭，定义了内积的线性空间称为内积空间
当K为数域时，（1）为(u,v)=(v,u)；
若(u,v)=0，则称u与v正交

柯西-施瓦茨不等式：

|(u,v)|²≤(u,u)(v,v)

格拉姆矩阵：

设X是一个内积空间，μ₁,μ₂,···,μ_n∈X，矩阵

称为格拉姆矩阵，G非奇异的充分必要条件是μ₁,μ₂,···,μ_n线性无关
该内积空间可以导出一种范数||u||=√((u,u))

权函数：

设[a,b]是有限或无限区间，非负函数ρ(x)满足：
（1）∫_a^bx^kρ(x)dx存在且为有限值
（2）当且仅当g(x)恒等于0时，∫_a^bg(x)ρ(x)dx=0
此时ρ(x)称为权函数

3.1.4 最佳逼近

在函数逼近问题中，我们把逼近的误差记为max||f(x)-P(x)||，我们一个让我们的逼近最有效的方法就是让max||f(x)-P(x)||最小化，也就是min max||f(x)-P(x)||，这时所得到的P(x)称为最佳逼近多项式。当我们使用范数去衡量误差大小时，当范数取|| ||_∞时，P(x)称为最优一致逼近多项式，当范数取|| || ₂时，P(x)称为最佳平方逼近多项式。