第一章 基本运算

上必有一根。
y f (x)
a x*
b
x
二分法的算法实现
a
xa1
x* xb2 b
给定有根区间 [a, b] ( f(a) f(b) < 0) 和 精度 1. 令 x = (a+b)/2 2. 如果 b – a < ， 停机，输出 x
3. 如果 f (a) f (x) < 0 , 则令 b = x，否则令 a = x, 返回 第1步
例1： 利用牛顿迭代法求解 f(x)=ex-1.5-tan-1x 的零点。初始 点 x0=-7.0
例1： 利用牛顿迭代法求解 f(x)=ex-1.5-tan-1x 的零点。初始
点 x0=-7.0
解： f (x0)=-0.702×10-1，f '(x)=ex-(1+x2)-1
计算迭代格式:
xk 1
给定一组离散数据： x x1 …… xn y y1 …… yn
要求一个函数 p(x) ，按照最小二乘原理，使得
达到最小值，这一方法被称为拟合，其中 p(x) 称为拟合 曲线。
一种常见的拟合曲线为多项式拟合，即选取拟合曲线 p(x) 为m次多项式
则问题变为求 ak (k=0,1,2,…,m), 使得 极小。
xk
f ( xk ) f ( xk )
它们的原函数都不是初等函数
数值积分的出发点
积分转化为求和
将区间 [a,b] 分割为 n 等份，每个小区间 的宽度为 h=(b-a)/n
线性近似——梯形法则
在每个区间 [xi, xi+1] 进行线性插值
二阶多项式近似——辛普生法则
在区间 [xi-1, xi+1] 上对 f(x) 进 行二阶插值。
注意
➢ 不同公式的精度：五点>三点>两点 ➢ 注意误差随步长 h 的减小先减小再增大
舍入误差
微分公式涉及两个很接近的 f 值相减。当步长过小时， 计算机的舍入误差会使导数计算不准确。
高阶导数
二阶导数
这也容易从二阶导数的定义直接得出 更高阶的阶导数以此类推
1.4 数值积分
为什么要数值积分？
牛顿-莱布尼兹公式
几何意义
y
x*
x
x0
x1
x0
f ( x0 ) f ( x0 )
迭代形式为
xk1 xk
f ( xk ) f ( xk )
牛顿法的算法构造
xi1
xi
f ( xi ) f ( xi )
[牛顿迭代法] 1: 初始化 x0 , δ,置 i:=0 2: 如果| f(xi ) |≤δ ，则停止. 3: 计算 xi+1:=xi - f (xi) / f ' (xi) 4: 如果 | f (xi+1) |≤δ ，则停止. 5: i:=i+1, 转至3.
不同积分方法的结果比较
注意，此时误差随h减小而减小，舍入误差并不重要，这是 因为积分公式中，所有 f 的值的符号都一样。
高阶的算法
选取更多的点，进行更高阶的插值可以得到更高 阶的算法，如
simpson 3/8 算法（三阶插值）
Bode规则（四阶插值）
过高阶的插值可能导致严重的振荡行为，从而给出被积 函数不准确的插值。所以为了得到更高精度，往往用低 阶插值，同时减小 h。
已知离散数据
x x0 x1 …… xn y y0 y1 …… yn
n
插值多项式为 y(x) Aj (x) y j j0
其中 满足
Aj (x)
n i0
x xi x j xi
i j
这称之为拉格朗日多项式插值。
讨论
是否阶数越高，效果越好？
例：连续函数
f
(
x)
1
5 x
2
L10 (x)
f (x)
在区间[-5,5]上取等距插值节点
-5
5
可以看出，L10(x)的误差在区间两端非常大
过高阶的插值可能导致严重的振荡行为，即Runge现象。
怎样改进？
解决的办法——分段插值
用分段低次插值，最简单的就是分段线性插值
不光滑！
样条插值
解决之道——在每个区间上，不是用线性函数而是用三 阶多项式进行插值。
区间[a,b] 有离散点：a= x0< x1< …< xn=b，及其对应的函 数值 yi ( i=0,1,…,n) ，如果函数S(x)满足条件：
(1) 在[a,b]上具有二阶连续导数，即 (2) 在每个子区间[xi , xi+1] ( i=0,1,…,n-1)上是三次多项式
(3)
则称S(x)是y = f(x) 的三次样条插值函数
确定S(x)需要4n个条件, 而 我们只给出了4n-2个条件
端点函数值：2个
内节点函数值及连续条件：2(n-1) 个 内节点一、二阶导数连续条件：2(n-1) 个
(x0,y0)
x0
x1
X
xn
寻找一个解析形式的函数 φ(x)， 满足 φ(xi) = yi, i=1,2,…, n
最常用的函数形式是多项式，称为多项式插值。
最直观的方法——以二阶为例
离散数据为
x y
x0 x1 x2 y0 y1 y2
设插值函数为 2 (x) a0 a1x a2 x2
带入数据，得
二分法的优缺点
✓ 简单易用 ✓ 稳妥，对 f (x) 要求不高，只要连续即可收敛
✓ 收敛速度慢
二分法——例题
例1: 用二分法求方程 2x3 5x 1 0 在区间 (1,2)内
的实根, 要求误差限为 102 。
5
4
f (x) 2x3 5x 1
3
2
1
0
-1
-2
-3
-4 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2
高于四次方程的一般代数方程没有 一般形式的代数解
f ( x) an x n an1x n1 a1x a0 (an 0)
更不用说更为复杂的方程
cos(x)cosh(x) 1 0
阿贝尔(1802—1829)
1.5.1 二分法（搜索法）
求 f (x) = 0 的根
原理：若 f C[a, b]，且 f (a) · f (b) < 0，则 f (x)=0 在 (a, b)
第一章 基本数学运算
本章要学会计算的物理问题：双原子分子的振动能级
本章内容
1 2 3 4 54 64
插值 拟合 数值微分 数值积分 方程求根 分子振动的半经典量子化
1.1 插值
问题的提出
给定一组离散的数据 x x0 x1 x2 …… xn y y0 y1 y2 …… yn
Y
(x1,y1)
(xn,yn)
3 -0.5
4
5
-0.25 0
6
7wenku.baidu.com
0.25
0.5
8
9
0.75
1
yi
-0.2209 0.3295 0.8826 1.4392 2.0003 2.5645 3.1334 3.7601 4.2836
设所求的最小二乘二次拟合多项式为
相应的正则方程为
其解为a0 = 2:0034; a1 = 2:2625; a2 = 0:0378，所求多 项式为
a
I b f ( x)dx F (b) F (a)
但是这要求
➢被积函数 f(x) 有解析表达式 ➢f(x) 的原函数 F(x) 为初等函数
我们面临的的问题
1) f(x) 没有解析表达式
x 0.1 0.2 0.3 0.4 0.5 f(x) 4 4.5 6 8 8.5
2) f(x)有解析表达式，但原函数不是初等函数 ，例如
由多元函数取极值条件，得 这就是正则方程。解这个方程即可得到拟合函数的形式。
当最高幂次为1时，即为最小二乘法的拟合直线，正 则方程简化为
Matlab 指令
polyfit(X, Y, N)
X, Y N 返回值
数据点及其函数值 拟合多项式的最高幂次 为按降幂排列的多项次系数
例子
i
1
xi
-1
2 -0.75
a0 a0
a1 x0 a1 x1
a2 x02 a2 x12
y0 y1
a0 a1 x2 a2 x22 y2
解这个方程，即可得到相应系数 a0, a1, a2。
N阶呢
离散数据为 x x0 x1 …… xn y y0 y1 …… yn
设插值函数为n阶多项式 n (x) a0 a1x a2 x2 an xn
1.3 数值微分
为什么要学习数值微分？
我们碰到的函数往往没有解析形式，例如可能是如下数表形式
x 0.1 0.2 0.3 0.4 0.5 f(x) 4 4.5 6 8 8.5
这就需要借助于数值微分 更重要的，数值微分是其它很多数值方法的基础
已知
x … -2h -h 0 h 2h …
f(x) … f-2 f-1
ti 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Ii 3.16 2.38 1.75 1.34 1.00 0.74 0.56
关于拟合，要提醒的
With four parameters I can fit an elephant, and with five I can make him wiggle his trunk.
带入数据，得
a0 a0
a1 x0 a1 x1
a2 x02 a2 x12
an x0n an x1n
y0 y1
a0 a1xn a2 xn2 an xnn yn
当n很大时，求解这个方程计算量太大，需另寻它法
拉格朗日插值——从一阶说起
对离散数据 x x0 x1 y y0 y1
插值方程为一直线方程，可表示为
f (x1) f (1.5)<0 得 …….
I2=[1.5, 1.75] , x2 =(1.5+1.75)/2=1.625
x*x7 =1.6758
I6=[1.681875, 1.6875], I7=[1.671875, 1.679688] b7 - a7=0.781310-2 < 10-2
1.5.2 牛顿法
f0
f1
f2
…
求 X=0 处的导数
方法：泰勒展开
x=h 处的函数值 x=-h 处的函数值
两点公式
向前差分 向后差分
三点公式（中心差分公式）
进一步的改进——5点公式
虽然精度更大，但是计算量也更大！ 通常的计算三点公式已足够。
例子：不同微分方法的精度比较
不同方法计算 d sin x/dx |x=1 = 0.540302 的误差
二分法——例题
例1: 用二分法求方程 2x3 5x 1 0 在区间 (1,2)内
的实根, 要求误差限为 102 。
解：令 f (x) 2x3 5x 1
f (1)<0, f (2)>0 记 I0=[1,2] , x0 =(1+2)/2=1.5
因为 f (x0) f (1)>0 得 I1=[1.5, 2] , x1 =(1.5+2)/2=1.75
反常积分的处理
反常积分分为两类：
➢ 积分区间有限,在积分区间内被积函数有奇点 ➢ 积分区间为无限
策略 通过积分变量的变换，将反常积分变换为普通积分
（1） 积分区间内含有奇点的积分
在 x=1 处有一个奇点，假设函数 g 在这点的值有 限，则积分为有限值。 做变换 t=(1-x)1/2，积分变为
（2）无限区间的积分
为了推广到高阶，将其写成更对称的形式
其中
A0
x x1 x0 x1
A1
x x0 x1 x0
满足
函数
节点 x0
x1
A0 (x)
10
A1 (x)
01
更进一步——二阶插值
y(x) A0 (x) y0 A1(x) y1 A2 (x) y2
Y y1
函数
节点 x0
x1 x2
yy02
A0 (x)
1 00
需补充2个边界条件，根据实际情况有不同的取法 例如，可取周期边界条件：
样条插值的例子
Matlab 指令
interp1(X, Y, Xi， method)
X, Y Xi method 返回值
数据点及其函数值 待求的数据点 插值方法，可取’linear’,’spline’等 数据点Xi的函数值
1.2 拟合
g(x) 在 x 很大时趋于常数，积分为有限值。 做变换 t=x-1，积分变为
Matlab自带的积分指令
➢ quad 用自适应辛普森算法。根据积分精度的需要， 自动调节积分取点的数目。
➢ 调用格式为 quad(函数句柄, 积分下限，积分上限)
quad(@(x)sin(x), 0.5, 0.6)
1.5 数值求根
y=f(x)
A1 (x) A2 (x)
0 10
0
01
O
x0
X
x1
x2
A0
(x ( x0
x1)(x x2 ) x1)(x0 x2 )
A1
(x (x1
x2 )(x x0 ) x2 )(x1 x0 )
A2
(x ( x2
x0 )(x x1) x0 )(x2 x1)
一般的——N阶多项式插值
实例——密立根油滴实验
指数拟合
如果数据点的分布近似于指数曲线，则可考虑用指数函数
去拟合数据。 对上式取对数
同时对数据 yi 也取对数得 Inyi, 利用数据组 (xi, Inyi) 求出 最小二乘拟合直线 。再取指数即得到原数据组的最小拟 合指数
例子
已知一发射源的发射强度具有指数形式 I=I0e-αt, 现有 一组观测数据如下