科学计算的背景与概况
数值分析是研究科学计算中各种数学问题求解的数值计算方法法
数值分析、科学计算与数值计算是当今科学研究的三种基本手段之一,是数学和计算机应用于其他学科的桥梁,本课程只考虑软件方面即数值计算的有关算法(数值仿真软件),数值分析和数值计算是研究求解连续数学问题(对象)的算法(核心)的学科
数值分析
概念:数值分析是研究各种数学问题的数值方法的设计、分析、有关的数学理论和具体实现的一门学科,实际上就是介绍用计算机解决数学问题的计算方法及其理论,利用计算机高速的简单运算去实现各种复杂的功能
地位:现代科学研究的三大支柱是科学理论、科学实验和科学计算,而科学计算的核心内容就是以现代化的计算机及数学软件为工具,以数学模型为基础进行模拟研究
在建立了数学模型之后,并不能立刻用计算机直接求解,还必须寻找用计算机计算这些数学模型的数值方法,即将数学模型中的连续变量离散化,转换成一系列相应的算法步骤,编制出正确的计算程序,再上机计算得到满意的数值结果
特点:
- 与计算机应用密切结合的实用性很强的学科
- 既要讨论 变量问题又要讨论离散变量问题,关心的是数值结果
- 近代数学的一个重要分支,专门研究数学问题的数值解法
研究内容:
- 插值问题
- 函数逼近
- 数值积分与数值微分
- 线性代数方程组的数值解法
- 非线性方程组的数值解法
- 代数特征值问题
- 常微分方程的数值解法
研究方法:递归迭代、近似代替、离散化、外推法
评判标准:误差、稳定性、收敛性、计算量、存储量和自适应性
科学计算与数值算法
数值计算的特点:处理连续数学的量(实数量),问题中常涉及微分、积分和非线性。被求解的问题一般没有解析解、或理论上无法通过有限步计算求解,分为无解析解、有解析解等问题,目标是寻找快速结束的算法
好数值算法的特点:
- 计算效率高、计算复杂度低
- 可靠性好:在考虑实际计算的各种误差情况下,结果尽可能的准确
数值计算的步骤:
- 建立数学模型
- 研究数值求解方程的算法
- 通过计算机软件实现算法
- 在计算机上运行软件进行数值模拟
- 将计算结果用较直观的方式输出,如图形可视化方法
- 解释和验证计算结果,如果需要重复上面的某些步骤
设计数值方法的关键:将问题简化(估计带来的误差),然后求解简化后的问题
数值计算中的误差
来源与分类
- 从实际问题中抽象出数学模型是产生的误差——模型误差
- 通过测量得到模型中参数的值,导致输入数据的误差——观测误差
- 求解近似——截断误差(方法误差)
- 机器字长限制——舍入误差
由于前两种误差不可避免,因此在科学计算中主要研究截断误差和舍入误差对计算结果的影响
绝对误差、相对误差
绝对误差:近似值 x 与准确值 x* 的差值 e
,误差一般不能准确计算,只能根据测量或计算情况估计出它的绝对值的一个上限,这个上限成为近似值 x 的绝对误差限,简称误差限,记为 ε,即
其意义是 -ε ≤ x-x* ≤ ε。在工程中常记为 x* = x ± ε
厘米刻度尺去量物体的时候,如果边缘处在两个刻度值的中间,那么测量值和真实值也仅仅最多相差 0.5cm,所以误差限就是 0.5cm
注意:绝对误差有时并不能很好的反应近似程度的好坏,10000内差 5 显然比 10 内差 1 要好,所以生活中常用相对误差来作为评判标准
相对误差:理论上是绝对误差与准确值 x* 的比值,即
相应的,若正数 εr 满足则如下关系式称 εr 为 x 的相对误差限
相对误差和绝对误差之间如下满足关系式
实际上由于 x* 不知道所以常用 x 代替 x* 作为分母,后来也就是用如下关系式计算相对误差了
有效数字
注意数值分析中的有效数字和高中学过的有效数字的定义是不同的,可以看这篇 参考博客
定义:如果 x 的误差限 ε 是它某一数位的半个单位,我们就说 x 准确到该位,x 从这一位起直到前面第一个非零数字为止的所有数字称为 x 的有效数字,即
则说 x 近似表示 x* 准确到小数后第 n 位,并且 x 从这小数点后第 n 位起直到最左边的非零数字之间的一切数字称为有效数字,并把有效数字的位数称为有效数位
注意说法,“近似表示”说明小数后这 n 位不一定相同,而是两者差的绝对值在误差范围内。无疑近似值的有效数字越多,其计算结果的绝对误差就越小
例如
π=3.1415926535...
,如果近似值为 3.14 则 |e|= |π-3.14|=0.0015926,小于等于 0.5*10-2,所以精确到小数点后两位有效数字就是 3.14,有效数位为 3。同理如果近似值为 3.1416 则 |e|=|π-3.1416|=0.0000074,小于等于 0.5*10-4,所以精确到小数点后四位有效数字就是 3.1416,有效数字为 5。关键是要会根据 |e| 和 准确值确定近似值的有效为主
定理:设近似值 x 有 n 为有效数字,则其相对误差限满足如下关系式,并且至少有 n 位有效数字
注意问题:计算 1/795 - 1/760
,视已知数为确切之,用四位浮点数计算
- 如果将两个分数分别表示为 4 位有效数字的浮点数直接计算,最终答案得 0.2*10-5,只剩下一位后效数字,损失大量有效数字造成相对误差的扩大
- 如果先将两个分数通分表示为 4 位有效数字的浮点数再计算,最终答案得 0.1734*10-5,得到 4 位有效数字与上面方案相比更加合理
造成这个的原因就是误差的传播与积累,生活中常见的就是蝴蝶效应的病态问题
数值计算原则
要使用数值稳定的算法
要避免两个相似数的相减
在数值计算中,两个相近的数作减法是有效数字会损失
对于两个相近的数相减,若是找不到适当的方法代替只能在计算机上采用双倍字长计算来提高精度
绝对值太小的数不宜作除法
-
由于除数很小,将导致商很大,很可能出现数据溢出
-
另外设 x* 和 y* 的近似值为 x,y,则
z=x/y
是 z* 的近似值。此时 z 的绝对误差满足如下公式,可见除数太小会导致商的绝对误差很大
例如
2.1782/0.001=2718.2
,即使分母变化很小为 0.0011,也会导致结果变为2.1782/0.0011=2471.1
避免大数吃小数
求和计算时避免运算符左侧放大数,右侧放小数。选择从小到大相加可以使和的误差减小
正确的方法是求出 x1
后,为了减少大数+小数造成的误差借助韦达定理求 x2
先化简再计算
目的是为了减少近似值的计算步骤从而避免误差的过度积累
课后作业
1 |
|
1 |
|
差异:公式一计算的数值非常不稳定,尤其是 I13 及以后的值几乎脱离正确答案,而公式二计算的数值更加稳定,一直到 I0 答案也贴近准确值
分析:
公式一中第 n 项由 1 减去乘上整数 n 的前一项结果递推而来,每次计算都会将前面各项产生的舍入误差扩大 n 倍并累加到当前项,另外由于项数逐渐增大每次计算的误差也放得越来越大,最终不断递归造成后续项结果严重失真
而公式二中第 n-1 项由 1 减去后一项的结果并除 n 递推而来,每次计算都会将前面各项产生的舍入误差缩小 n 倍,因此即使递归多次误差也不会明显的扩大从而保证了每一项都接近准确值
非线性方程求根
对于方程 f(x)=0
,若 f(x) 是一次多项式则称为线性方程,否则成为非线性方程
非线性方程分为两种:
- 代数方程:a0xn + a1xn-1 + … + an-1x + an = 0,其中 a0 ≠ 0,ai ∈ R,例如 x3 - x - 1 = 0
- 超越方程:例如 x-e-x = 0
对于代数方程当 n=1,2,3,4
可以使用相应的求根公式,但是对于 n>=5
则不存在求根公式了
非线性方程一般可能有无穷多个解,求解时需要强调求解区间,非线性方程一般没有直接解法,通常二分法确定大致区间,使用迭代算法求解
方程
f(x)=0
的根称为函数 f(x) 的零点。如果 f(x) 表示为 f(x) = (x - )mg(x),g(a) ≠ 0,则称 a 为f(x)=0
的重根。当 m = 1 时称为单根,当 m > 1 时称为方程的 m 重根,此时有 f(a) = f’(a)= … = f(m-1)(a) = 0,fm(a) ≠ 0
二分法(对分法)
理论介绍
基本思想:将有根区间对分,并找出根所在的小区间,然后再对小区间对分直到有根区间的长度足够小为止
数学原理:零点定理——设 f(x) 在 [a,b]
上连续并且 f(a)f(b)<0
则在区间 (a,b)
内至少存在一点 a 使得 f(a)=0
适用范围:求有根区间内的单根或者奇重实根,即 f(a)f(b)<0
成立的区间
算法步骤:
- 计算 f(a) 和 f(b),若
f(a)f(b)>0
则算法失效停止计算 - 另 x 取中间的值
x=(a+b)>>1
,计算 f(x) - 若是
|f(a)|<eps
或者|b-a|<eps
则停止计算输出近似解 x - 若
f(a)f(x)<0
则令b=x
否则令a=x
- 返回第二步
eps 为可接受的最大误差范围
优点 | 缺点 |
---|---|
简单易用,总是收敛 | 收敛慢,不能永安里求偶数重根,一次只能求一个根 |
总结:一般用来计算解的一个粗糙估计
收敛性证明与误差分析:
结论:二分总是收敛的并且误差估计式为
因为不知道准确值 x* 所以误差限不能用前面的方法得到,但是根据导出的误差估计式可以知道
代码模板
1 |
|
例题
求 x3 - x - 1 = 0 在 [1.0,1.5]
内的一个实根,准确到小数点后两位
根据误差估计式,想要精确到小数点后两位,则需要满足 |x* - x6| ≤ 0.05,也就是 |b-a| / 2k-1 ≤ 0.05,根据表可以看出需要 6 次二分
在实用中常用二分法来判别根的存在区间,如区间较大可用二分法适当收缩区间,并选择初值 x0 为该区间中点,再用收敛速度快的迭代法迭代计算求根
不动点迭代法
理论介绍
基本原理:将方程 f(x)=0
改写为等价形式 x=Φ(x)
。若是 x* 满足 f(x*)=0
,则 x*=Φ(x*)
,反之依然我们称 x* 为 Φ(x) 的一个不动点。求 f(x) 的零点就等价于求 Φ(x) 的不动点,选择一个初始近似值 x0 作为起点,将它带入 x=Φ(x)
的右端即可求出 x1,如此反复迭代计算就可以求出 xk+1 = Φ(xk) (k=0,1,2,…),Φ(x)
被称作迭代函数
如果对任何 x0 ∈ [a,b] 得到的序列 {xk} 有极限
则称迭代方程收敛,且 x* = Φ(x*) 为 Φ(x) 的不动点,称此通过收敛迭代函数求 x* 的方法为不动点迭代法
数学原理:上述迭代法是一种逐次逼近法,其基本思想是将隐式方程 x=Φ(x)
归结为一组显式的计算公式 xk+1 = Φ(xk),实质上就是一个逐步显示化的过程,方程 x=Φ(x)
的求根问题在 xOy 平面上就是要确定曲线 y=Φ(x)
与直线 y=x
的交点 P*,对于 x* 的某个近似值 x0,在曲线 y=Φ(x)
上可以确定一个点 P0,它以 x0 为横坐标,而纵坐标则等于 Φ(x0)=x1,过 P0 引平行 x 轴的直线,设此直线交直线 y=x
于 Q1 点,然后过 Q1 再作平行于 y 轴的直线,与曲线 y=Φ(x)
的交点记作 P1,则点 P1 的横坐标为 x1,纵坐标为 Φ(x1)=x2
如上图箭头所示的路径依次做下去求得迭代值 x1,x2,…,如果点列 {Pk} 趋向于 P* 则相应的迭代值 xk 收敛到所求的根 x*
注意:迭代函数可能不止一个,变式方法不同得到的也不同,其中有的收敛有的发散
例如求方程 f(x) = x3 - x - 1 = 0 在 x0 = 1.5 附近的根 x*,如下两种迭代函数第一种收敛,第二种发散
如果仅取 6 位数值,对于第一种收敛的迭代函数我们迭代 7 次就可以找到所要的根
但是第二种发散的迭代函数结果会越放越大不可能趋向于某个极限值,纵使迭代千百次结果也是毫无意义的,所以我们需要的迭代函数是收敛去 P* 的
不动点的存在性与迭代法的收敛性:
局部收敛与收敛阶:
模板代码
1 |
|
牛顿法
理论介绍
基本思想:一种线性化思想,将非线性方程 f(x)=0
逐步归结为某种线性方程来求解
迭代公式:
几何解释:方程 f(x)=0
的根 x* 可以解释为曲线 y=f(x)
与 x 轴的交点的横坐标,设 xk 是根 x* 的某个近似值,过曲线 y=f(x)
上的横坐标为 xk+1 作为 x* 的新的近似值,由于这种几何背景牛顿法亦称为切线法
收敛性:牛顿迭代法的主要优点是收敛速度比不动点迭代法还要快,为平方收敛
牛顿法的两个缺点,其中一个是每次迭代都要重新计算导数,计算量太大,另一个就是收敛性依赖于 x0 的选择
计算步骤:
模板代码
1 |
|
牛顿法收敛块但是确实也很明显,一是每步迭代都要计算 f(xk) 和 f’(xk) 计算量较大且有时计算导数过于困难,二是初始近似值 x0 值在根附近是才能保证收敛,为了克服上述缺点通常采用下述改进的牛顿法
改进牛顿法
简化牛顿法
理论介绍
基本思想:用 f’(x0) 直接替换所有的 f’(xk),简化牛顿法也称为平行弦法
迭代公式:
该迭代法在根 x* 的附近局部线性收敛,因此比原先的牛顿法收敛速度要慢,几何意义是用斜率为 f’(x0) 的平行弦与 x 轴的交点作为 x* 的近似
改进:解决了牛顿法每步迭代都需要计算导数的问题
模板代码
1 |
|
牛顿下山法
理论介绍
基本思想:牛顿法的 x0 如果偏离所求根 x* 较远那么会导致发散,要求每步迭代都满足下降条件,所以附加一项要求保证全局收敛的条件,即函数单调性满足:|f(xk+1)| < |f(xk)|
迭代公式:结合牛顿法的公式但是每次求出 xk+1 后与前一步的近似值 xk 适当加权平均改进当前的 xk+1
即采用如下公式求出改进的 xk+1,如果不满足单调性重新调整 λ(下山因子)
λ 的取法:从 1 开始逐次减半,直到满足下降条件为止
改进:解决了牛顿法迭代过程中发散的问题
模板代码
1 |
|
为了避开遇见导函数比较困难的情况还可以利用已求的 f(xk) 和 f(xk-) 直接回避计算导数,下述的弦截法和抛物线法就是此方法
弦截法
理论介绍
基本思想:设 xk 和 xk+1 是 f(x)=0
的近似根,利用 f(xk) 和 f(xk-1) 构造一个插值多项式 p1(x),并用 p1(x) = 0 的根作为新的近似根 xk+1
迭代公式:
弦截法和简化牛顿法(切线法)都是线性化方法,但是两者存在本质的区别,切线法计算 xk+1 只用到前一步的值 xk 但是弦截法需要用到前两步的值 xk 和 xk-1 所以这个方法必须起初选定两个初始值
模板代码
1 |
|
抛物线法
基本思想:用二次曲线与 x 轴的交点作为 x* 的近似值
收敛性:在一定的条件下可以证明抛物线的收敛阶为 p 约等于 1.840,抛物线是超线性收敛,收敛的速度比弦截法更接近牛顿法
课后作业
题目:用二分法求方程 x2 - x -1 =0 的正根,要求误差小于 0.05
分析:相当于近似精确到小数点后 1 位,手算一下大概处于区间 [1,2]
之间,根据二分法的误差估计式终止条件需要 1 / 2k+1 <= 0.05,最终计算结果为 1.59375,经验证与 (√5+1)/2=1.6180 处于误差范围之内
1 |
|
题目:比较求 ex + 10x - 2 = 0 的根保留三位有效数字所需的计算量:
- 在区间
[0,1]
借助二分法 - 用迭代法 xk+1 = (2 - ek) / 10,取初值 x0 = 0
分析:设置误差为 0.5*10-4 从小数点后第四位向前得到三位有效数字。
- 二分法根据误差估计式:终止条件需要满足 1 / 2k+1 <= 0.5*10-4
- 迭代法根据误差估计式:L = |Φ’(0)| <= 0.825,终止条件需要满足 L * |xk+1 - xk| / (1 - L) <= 0.5*10-4
计算量以中间计算的 x 个数为基准,经过计算统计二分法需要计算 15 次 x 而迭代法只需要计算 6 次 x,最终计算结果保留三位有效数字均是 0.0905 经验证正确
1 |
|
1 |
|