第一章 基本运算(计算物理matlab)

A1
( x x2 )( x x0 ) ( x1 x2 )( x1 x0 )
A2
( x x0 )( x x1 ) ( x2 x0 )( x2 x1 )
一般的——N阶多项式插值
已知离散数据
x y
x0 x1 y0 y1
n
……
……
xn yn
插值多项式为 y( x) Aj ( x) y j
解决的办法——分段插值
用分段低次插值，最简单的就是分段线性插值
不光滑！
样条插值
解决之道——在每个区间上，不是用线性函数而是用三 阶多项式进行插值。
区间[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)
j 0
其中 满足
x xi A j ( x) i 0 x j xi
n i j
这称之为拉格朗日多项式插值。
讨论
是否阶数越高，效果越好？ 例：连续函数
f ( x) 5 1 x2
L10 ( x)
f ( x)
在区间[-5,5]上取等距插值节点
-5
5
可以看出，L10(x)的误差在区间两端非常大 过高阶的插值可能导致严重的振荡行为，即RungeBiblioteka Baidu象。 怎样改进？
(xn,yn)
xn yn
X
xn
x0
x1
寻找一个解析形式的函数 φ (x)， 满足 φ (xi) = yi, i=1,2,…, n 最常用的函数形式是多项式，称为多项式插值。
最直观的方法——以二阶为例
离散数据为 x x0 x1 x2
y
设插值函数为 带入数据，得
y0 y 1 y2
2 ( x) a0 a1 x a2 x 2
反常积分的处理
反常积分分为两类：
积分区间有限,在积分区间内被积函数有奇点 积分区间为无限
策略
通过积分变量的变换，将反常积分变换为普通积分
（1） 积分区间内含有奇点的积分
在 x=1 处有一个奇点，假设函数 g 在这点的值有
限，则积分为有限值。 做变换 t=(1-x)1/2，积分变为
（2）无限区间的积分
牛顿法的算法构造
f ( xi ) xi 1 xi f ( xi )
[牛顿迭代法]
1: 初始化 x0 , δ,置 i:=0
求
X=0 处的导数
方法：泰勒展开
x=h 处的函数值
x=-h 处的函数值
两点公式
向前差分
向后差分
三点公式（中心差分公式）
进一步的改进——5点公式
虽然精度更大，但是计算量也更大！
通常的计算三点公式已足够。
例子：不同微分方法的精度比较
不同方法计算 d sin x/dx |x=1 = 0.540302 的误差
2. 如果 b – a < ， 停机，输出 x 3. 如果 f (a) f (x) < 0 , 则令 b = x，否则令 a = x, 返回
第1步
二分法的优缺点
简单易用
稳妥，对 f (x) 要求不高，只要连续即可收敛
收敛速度慢
二分法——例题
例1: 用二分法求方程 2 x 3 5x 1 0 在区间 (1,2)内
I

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