《计算物理学-李录》由会员分享，可在线阅读，更多相关《计算物理学-李录（170页珍藏版）》请在金锄头文库上搜索。
1、 计算物理学 李 录 山西大学物理电子工程学院 二 00 六年二月 目 录 i 目 录 第一章 基本数学运算 1 第一节 插值 1 第二节 拟合 3 第三节 数值微分 7 第四节 数值积分 9 第五节 求根 15 第二章 常微分方程的初始问题 23 第一节 几种简单的数值解法 23 第二节 Runge Kutta 方法 28 第三节 多步法 36 第四节 稳定性 40 第五节 动力学中的有序和混沌 44 第三章 边值问题和本征值问题 51 第一节 物理学中出现的边值问题和本征值问题举例 51 第二节 Numerov 算法 52 第三节 边值问题的直接积分 54 第四节 边值问题的 Green
2、函数法 57 第五节 打靶法 64 第六节 一维 Schr dinger 方程的定态解 69 第四章 特殊函数和 Gauss 求积法 74 第一节 特殊函数 74 第二节 Gauss 求积法 81 第三节 量子散射的 Born 近似和程函近似 86 第五章 矩阵中的数值计算方法 96 第一节 物理学中的矩阵 96 第二节 矩阵的基本运算 98 第三节 一般矩阵的本征值问题 104 第四节 对称矩阵的本征值问题 108 第六章 椭圆型微分方程 117 第一节 离散化和变分原理 117 第二节 求解边值问题的一种迭代方法 119 第三节 关于离散化的进一步讨论 127 第七章 抛物型微分方程 13
3、2 第一节 简单的离散化和条件稳定性 132 第二节 隐式格式和回代法 136 第三节 高维扩散 144 第四节 本征值问题的迭代方法 148 第五节 含时间的Schr dinger方程 155 目 录 ii 第八章 Monte Carlo方法 158 第一节 Monte Carlo 方法的基本思想 158 第二节 具有特定分布的随机变量得产生 162 第三节 Metropolis等人的算法 166 1 第一章 基本数学运算 本章将介绍数值计算中一些基本的运算 例如插值 拟合 数值微分 数值积分和求根等 第一节 插 值 Interpolation 当我们要从一组不完全的或离散的数据中获取某些局
4、部的信息时 需要使用函数的插值 也就是说 对于函数 yf x 通常是一个未知的或是一个较复杂的函数 xa b 若已知 a b 上一系列点 xxxxn 012 L处的函数值 yf x ii in 012 L 当我们要求这些点xxxxn 012 L之间的点x处的函数值时 就需用一个较简单的且满足条件 xy ii in 012 L 的函数 x来代替 f x 这一方法被称为插值法 其中 x被称为插值函数 f x被称为被插值函 数 点 i x被称为插值点 当这一较简单的函数 x被选择为多项式时 则称为代数插值 当这一较简单的函数 x被选择 为三角函数时 则称为三角插值 下面主要来介绍代数插值 1 线性插
5、值 也称一次插值 已知函数表 x x0 x1 y y0 y1 去构造一个线性多项式 1x 使其满足 1xy ii i 01 这样的一次插值多项式具有如形式 1 0 01 1 xy lxy l x 其中 lx xx xx 0 1 01 l x xx xx 1 0 10 被称为线性插值的基函数 误差为 R xf xxO h 1 2 其中 01 xxh 为步长 2 二次插值 已知函数表 x x0 x1 x2 y y0 y1 y2 去构造一个二次多项式 2x 使其满足 2xy ii i 012 这样的二次插值多项式为 xlyxlyxlyx 2211002 其中 lx xxxx xxxx 0 12 01
6、02 l x xxxx xxxx 1 02 1012 lx xxxx xxxx 2 01 2021 被称为二次 插值的基函数 误差为 R xf xxO h 2 3 其中h为步长 3 n次插值 2 已知函数表 x x0 x1 x2 xn y y0 y1 y2 yn 去构造一个n次多项式 nx 使其满足 n ii xy in 012 L n次插值多项式为 n i n i i xy l x 0 其中 l x xxxxxxxxxx xxxxxxxxxx i iin iiiiiiin 0111 0111 LL LL in 012 L 被称为n次插值的基 函数 误差为 R xf xxO h n n 1 其
7、中h为步长 下面是为n次插值的基函数而编写的 MatLat 函数文件 chap1 interpolyb m function yy interpolyb x1 x n yy 1 1 for i 2 n 1 yy 1 yy 1 x x1 i x1 1 x1 i end yy n 1 1 for i 1 n yy n 1 yy n 1 x x1 i x1 n 1 x1 i end for j 2 n yy j 1 for i 1 j 1 yy j yy j x x1 i x1 j x1 i end for i j 1 n 1 yy j yy j x x1 i x1 j x1 i end end 和
8、n次插值函数的函数文件 chap1 interpolyp m function y interpolyp x1 y1 x n a interpolyb x1 x n y a y1 说明当1 n和2 n时 分别可以得到一次插值函数和二次插值函数 例 在某实验中要测某短时光束的发光强度随时间的变化 已知在实验中测得的一组数据如下 t 为测数据的时间 y为反映光束强度的指标 求25 2 t秒处的光的强度指标 t 0 0 5 1 1 5 2 2 5 3 y 0 0 4794 0 8415 0 9975 0 9093 0 5985 0 1411 用六次插值多项式来求解该问题 运行的程序为 chap1 f
9、igure1 1 1 m x1 0 0 5 1 1 5 2 2 5 3 y1 0 0 4794 0 8415 0 9975 0 9093 0 5985 0 1411 3 n min length x1 length y1 n n 1 m 100 h max x1 min x1 m for i 1 m 1 x i min x1 i 1 h phin i interpolyp x1 y1 x i n end n interpolyp x1 y1 2 25 n plot x phin hold on plot x1 y1 bo 该程序的运行结果不仅可以显示出25 2 t秒处的光强度为 0 7781
10、同时可以画出光强度随时间的演 化 即插值函数 的图象 并描绘出它与原始数据之间的位置关系 如下图所示 00 511 522 53 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 1 第二节 拟 合 Fitting 正如上一节所讲 插值主要用于获取一个给定的离散集合的局部逼近 如果要获取一个给定的离散 集合的整体行为 我们则采用拟合 也就是说 给定一组离散数据 x 1 x 2 x 3 x n x y 1 y 2 y 3 y n y 要求一个函数 xpy 按照最小二乘原理 使得 2 1 ii n i yxp 达到最小值 这一方法被称为拟合法 其中 xp被称为拟合曲线 值
11、得注意的是拟合曲线 xpy 的确定可根据实验或实际情况来选择 通常我们考虑多项式拟合 即选取拟合曲线 xp为m次多项式 k k m k m xaxpxp 0 则问题变为求 2 1 0 mkakL 使得函数 图 1 1 1 光强度随时间的演化 4 2 01 10 i k ik m k n i m yxaaaafL 达到最小值 确定 2 1 0 mkakL 的方法是 让 mj a f j 2 1 0 0L 于是得到 mjxyax n i j iik n i jk i m k 2 1 0 010 L 这是一个1 m元线性方程组 特别地 当1 m时 即线性拟合 此时拟合直线 011 axaxp 由系数
12、 1 a和 0 a完全确定 它们 满足如下方程组 n i i n i i n i ii n i i n i i ynaxa yxxaxa 1 0 1 1 11 0 1 2 1 通过引入向量矩阵 1 1 1 2 1 n x x x A MM 和 n y y y b M 2 1 我们可以将该方程改写为如下的矩阵形式 bA a a AA TT 0 1 由此我们就可以得到拟合直线的系数为 bAAA a a TT 1 0 1 下面是为求拟合直线的系数 1 a和 0 a所编写的 MatLab 函数文件 chap1 fitlinec m function yy fitlinec x y nx length
13、x ny length y if nx ny warning The lengths of x and y should be same return end n min nx ny if n 2 error The number of the DATA should be greater than 1 return end x x 1 n y y 1 n x reshape x n 1 y reshape y n 1 A x ones n 1 5 b y B A A b A b yy B b yy yy 和函数文件 chap1 fitlinep m function p1 fitlinep
14、x y t a fitlinec x y p1 a t 1 一般地 我们引入矩阵 1 1 1 1 2 1 2 1 1 11 n m n m n m n m mm xxx xxx xxx A L MMOMM L L 和 n y y y b M 2 1 则拟合曲线 xpm的系数 011 aaaa mm L 所满足的方程组可以改写为 bA a a a AA Tm m T 0 1 M 由此我们就可以得到拟合曲线的系数为 bAAAaaa TT T mm 1 01 L 下面是为求拟合曲线 xpm的系数 011 aaaa mm L 所编写的 MatLab 函数文件 chap1 fitpolyc m func
15、tion yy fitpolyc x y m nx length x ny length y if nx ny warning The lengths of x and y should be same return end n min nx ny if n 2 error The number of the DATA should be greater than 1 return end if n restart l 0 t t x i x i 1 2 t x i 1 x i x i x i 1 2 x i x i 1 l 1 t t x i t x i 1 x i x i 1 2 x i
16、x i x i 1 2 x i 1 l 2 t t x i t x i x i 1 2 x i 1 x i x i 1 x i x i 1 2 simplify int l 0 t f i l 1 t f i 1 2 l 2 t f i 1 t x i x i 1 factor 1 6 4 f i 1 2 f i1 fi xix i1 这里 1122 1 1 ii xx iii xffffxff ii 下面我们举几个利用 MatLab 计算数值积分的例子 例 1 计算积分dxe x 1 0 显然 该积分的精确值为718281828 11 1 0 edxe x 下面我们分别利用梯形公式和 Simpson 公式来计算 利用梯形公式的 MatLab 命令文件程序 命令文件 chap1 trapezoid1 m 和函数文件 integrand1 m Using trapezoid rule to calculate the integral a 0 b 1 n 10000 h b a n int integrand1 a for i 2 n int int 2 integrand1 a i 1
