第 10 章计算物理与 Fortran 编程 - 习题课

第10章 计算物理与Fortran编程-习题课主讲教师：蒋臣威西安交通大学理学院

例题例题1. 利用数值方法求解f=sinx的一阶和二阶导数一阶导数的中心差分公式二阶导数的差分公式注意：对于边界点，可以采用向前差分、向后差分，或者外推法。（1）利用Matlab编程求解见Matlab程序 chap10_ex_1_sinx_differential.m （2）利用Fortran编程求解见Fortran程序 chap10_ex_1_sinx_differential.f90

数值积分——辛普生法则 在区间 [xi-1, xi+1] 上对 f(x) 进行Lagrange二阶插值。 2h

例题2：利用simpson方法计算 simpson法（1）利用Matlab编程求解见Matlab程序 chap10_ex_2_integration.m （2）利用Fortran编程求解见Fortran程序 chap10_ex_2_integration.f90 熟悉fotran程序的两种编译方式(a) ifort –c –r8 chap10_ex_2_integration.f90 ifort –o chap10_ex_2_integration chap10_ex_2_integration.o (更加精确) (b) Ifort –o chap10_ex_2_integration chap10_ex_2_integration.f90

3.1.方程数值求根的必要性 阿贝尔(1802—1829) 3. 方程求根阿贝尔-鲁芬尼定理：高于四次的一般代数方程没有一般形式的代数解对于更为复杂的方程，如非线性方程很难利用解析的方法求得方程的根。后面很多章节，如常微分方程的边值问题的打靶法求解等问题中也将用到方程的数值求根问题。作为基础我们只讨论一元非线性方程的数值求根问题。

3.2 二分法（对分搜索法） b x1 x2 x* a a b 求 f (x) = 0 的根根的存在性原理：若 fC[a, b]，且 f (a) · f (b) < 0，则 f (x)=0 在 (a, b) 上必有根。什么时候停止搜索？ y f (x) or a 如何搜索？ x* b x and

f(x) f(x) xk x* 只用可能造成误差偏大 xk+1 只用可能造成误差偏大 2 x* xk 二分法的算法实现给定有根区间 [a, b] ( f(a) f(b) < 0) 和精度要求1，2 1. 令 x = (a+b)/2; 2. 如果 |b – a| < 1并且|f(x)|< 2，结束运算，输出 x; 3. 如果 f (a) f (x) < 0 , 则令 b = x，否则令 a = x, 返回第1步

输入函数f(x),区间[a,b],误差标准1, 2 No if f(a)f(b)<0 Yes [a,b]上无根 if |a-b|> 1或 |f((a+b)/2)|> 2 No x=(a+b)/2 Yes x=(a+b)/2 n=n+1 end if f(x)f(a)>0 No Yes b=x a=x 二分法的流程图二分法的优缺点 • 简单易用 • 稳妥，对 f (x) 要求不高，只要连续即可收敛 X收敛速度慢，且只能求单根

3.2 Newton-Raphson方法 y x2 x x1 x* x0 设一元方程 f(x)=0的非线性函数f(x)连续可微。我们在其解的近似值xk附近将f(x)作Taylor展开几何意义因此可以构造迭代公式：什么时候停止迭代？

牛顿迭代法的算法构造 1: 初始化x0 , 误差标准δ,置 k=0 2: 如果| f(xk ) |≤δ，则停止. 3: 计算 xk+1=xk - f (xk) / f ' (xk) 4: 如果 | f (xk+1) |≤δ，则停止. 5: k=k+1, 转至3. 画出Newton算法的流程图

x0 x0 x0 算法说明注：Newton-Raphson Method 收敛性依赖于x0的选取。   x* • 收敛快，稳定性好，精度高等优点，是求解非线性方程的有效方法之一 X每次迭代均需计算函数值与导数值，计算量大。当导数值提供有困难时， Newton法无法进行。收敛性依赖初值选择。

3.4 弦割法 割线切线 x1 x0 弦割法在Newton-Raphon法的效率与必须计算导数的麻烦间提供了一种折中。向后差分切线斜率　　割线斜率代入Newton迭代公式任意2个初值 x0和 x1可以启动这个递推关系。

割线切线或 x1 x0 且 2 x* xk When to stop? x* 不能保证xk的精度

1: 初始化x0 、 x1, 误差标准δ与 ,置 k=0 2: 如果| f(xk ) |≤δ且| xk –xk+1 | ≤ ，则停止. 3: 计算 4: k=k+1, 转至2. 弦割法的算法构造画出弦割法的流程图

故，满足精度要求．在附近的根( ) 例3. 用弦割法求方程解：取由迭代公式求得下表见Matlab程序 chap10_ex_3_bisection_newton_secant.m （1）利用Matlab编程求解见Fortran程序 chap10_ex_3_bisection_newton_secant.f90 （2）利用Fortran编程求解

（1） 当关于和是线性时，式（1）为线性两点边值问题（2）常微分方程边值问题的数值解法在具体求解常微分方程时，必须附加某种定解条件。定解条件通常有两种，一种是初始条件，另一种是边界条件。与边界条件相应的定解问题称为边值问题。本节介绍求解两点边值问题的数值解法。

（3） （4）这里的为在处的斜率。令，上述二阶方程可降为一阶方程组打靶法打靶法的基本原理是将两点边值问题（1）转化为下列形式的初值问题

因此，边值问题变成求合适的 ，使上述方程组初值问题的解满足原边值问题的右端边界条件，从而得到边值问题的解。这样，把一个两点边值问题的数值解问题转化为一阶方程组初值问题的数值解问题。方程组初值问题的所有数值方法在这里都可以使用。问题的关键是如何去找合适的初始斜率的试探值。对给定的，设初值问题（3）的解为，它是的隐函数。假设随是连续变化的，记为，于是我们要找的就是方程求根函数F(x)

（5） 这样，可以按下面简单的计算过程进行求解。先给定两个初始斜率，分别作为初值问题（4）的初始条件。用一阶方程组的数值方法求解它们，分别得到区间右端点的函数的计算值和。如果或，则以或作为两点边值问题的解。否则用割线法（5）求，同理得到，再判断它是否满足精度要求。的根。可以用第2章的二分法或迭代法求上述方程的根。比如用弦割法有

如此重复，直到某个满足 ，此时得到的和就是边值问题的解函数值和它的一阶导数值。上述方程好比打靶，作为斜率为子弹的发射，为靶心，故称为打靶法。

打靶法的算法实现 1. 给定初始点的斜率猜测值S1； 2. 用常微分方程的初值问题解法（如RK算法）求解y(b,S1); 3. 给出另一个斜率猜测值S2作为弦割法的第二个启动点; 4. 用常微分方程的初值问题解法（如RK算法）求解y(b,S2); 5. IF (abs(y(b,S1)-y(b))<精度  ), y(x)=y(x,S1); 6. IF (abs(y(b,S2)-y(b))<精度  ), y(x)=y(x,S2); else 7. 利用弦割法迭代公式求出下一个斜率S3； 8. IF (abs(y(b,S3)-y(b))<精度  ), y(x)=y(x,S3); else S1=S2; S2=S3；回到第7步 9. 得到边值问题的解y(x)=y(x,S3);

例 4用打靶法求解非线性边值问题 要求误差不超过，其解析解是。解：对方程进行降阶处理，并转化成初值问题见Matlab程序 chapter10_ex_4_shooting_method.m （1）利用Matlab编程求解见Fortran程序 chapter10_ex_4_shooting_method.f90 （2）利用Fortran编程求解

随机数的生成-Metropolis算法 虽然前面我们讲述了多种按照规定分布产生随机数的有效方法，但要推广这些方法以对多维空间中的一个复杂的权函数抽样却是很困难或者是不可能的。一种很普遍的产生具有任意形状的给定概率分布的随机变量的方法叫做Metropolis -Rosenbluth-Rosenbluth-Teller-Teller 算法（简称Metropolis算法）。该算法已经广泛第应用在统计力学等很多学科中。 Metropolis (1915-1999) Nicolas Metropolis,Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, Edward TellerJ. Chem. Phys. 21, 1087 (1953). — The Paper was cited 14776 times from 1988 to 2013

Metropolis算法的实现 设想在变量X的空间(可能是多维的)内产生一组按概率密度 w (X )分布的点，Metropolis算法假想有一个随机行走者在X空间中运动，该随机行走过程相继各步的终点产生出点的一个序列：X0, X1,…；随着行走的路程越长，它连接的点就越接近所要求的分布。随机行走在位形空间中进行的规则为：设行走者处于序列中的Xn点上，为了产生Xn+1，行走者迈出试探性的一步到一新点 Xt,该新点可以用任何方便的方法选取。例如可以在点Xn周围的一个边长d很小的多维立方体中均匀地随机选取，然后按比值

来决定是“接受”还是“拒绝”该步试验。如果 r >1，接受该步(即取Xn+1=Xt)；而如果r <1，则以概率r 接受这一步。这时我们把 r和一个在［0,1］区间上均匀分布的随机数h比较，若h<r就接受这一步。如果该步试验未被接受，就舍弃，取Xn+1=Xn。这样产生出Xn+1之后，可以再从Xn+1出发迈出一个试验步按照同样的过程产生Xn+2。任意点X0都可用作随机行走的起点。上述算法的核心思想就是尽量多地选择w（X ）较大的点。

分布的样点 例题5用Metropolis算法产生按照见Matlab程序 chapter10_ex_5_metropolis.m （1）利用Matlab编程求解见Fortran程序 chapter10_ex_5_metropolis.f90 （2）利用Fortran编程求解

第 10 章计算物理与 Fortran 编程 - 习题课