计算物理学-李录

计算物理学

李 录

山西大学物理电子工程学院

二00六年二月

目 录

第一章 基本数学运算 .................................................................... 1

第一节 插值 ........................................................................ 1 第二节 拟合 ........................................................................ 3 第三节 数值微分 .................................................................... 7 第四节 数值积分 .................................................................... 9 第五节 求根 ....................................................................... 15 第二章 常微分方程的初始问题 ........................................................... 23

第一节 几种简单的数值解法 ......................................................... 23 第二节 Runge-Kutta方法 ............................................................ 28 第三节 多步法 ..................................................................... 36 第四节 稳定性 ..................................................................... 40 第五节 动力学中的有序和混沌 ....................................................... 44 第三章 边值问题和本征值问题 ........................................................... 51

第一节 物理学中出现的边值问题和本征值问题举例 ..................................... 51 第二节 Numerov算法 ............................................................... 52 第三节 边值问题的直接积分 ......................................................... 54 第四节 边值问题的Green函数法 ..................................................... 57 第五节 打靶法 ............................ ........................................ 64 第六节 一维Schrödinger方程的定态解 ................................................ 69 第四章 特殊函数和Gauss求积法 ......................................................... 74

第一节 特殊函数 ................................................................... 74 第二节 Gauss求积法 ............................................................... 81 第三节 量子散射的Born近似和程函近似 .............................................. 86 第五章 矩阵中的数值计算方法 ........................................................... 96

第一节 物理学中的矩阵 ............................................................. 96 第二节 矩阵的基本运算 ............................................................. 98 第三节 一般矩阵的本征值问题 ...................................................... 104 第四节 对称矩阵的本征值问题 ...................................................... 108 第六章 椭圆型微分方程 ................................................................ 117

第一节 离散化和变分原理 .......................................................... 117 第二节 求解边值问题的一种迭代方法 ................................................ 119 第三节 关于离散化的进一步讨论 .................................................... 127 第七章 抛物型微分方程 ................................................................ 132

第一节 简单的离散化和条件稳定性 .................................................. 132 第二节 隐式格式和回代法 .......................................................... 136 第三节 高维扩散 .................................................................. 144 第四节 本征值问题的迭代方法 ...................................................... 148

第五节 含时间的Schrödinger方程 ................................................... 155

第八章 Monte-Carlo方法 ..... .......................................................... 158

第一节 Monte-Carlo方法的基本思想 ................................................. 158 第二节 具有特定分布的随机变量得产生 .............................................. 162 第三节 Metropolis等人的算法 ....................................................... 166

第一章 基本数学运算

本章将介绍数值计算中一些基本的运算，例如插值、拟合、数值微分、数值积分和求根等．

第一节 插 值 (Interpolation)

当我们要从一组不完全的或离散的数据中获取某些局部的信息时，需要使用函数的插值．也就是说，对于函数y=f(x)（通常是一个未知的或是一个较复杂的函数），x∈[a,b]，若已知[a,b]上一系列点

x0,x1,x2,L,xn处的函数值

yi=f(xi)，i=0,1,2,L,n

当我们要求这些点x0,x1,x2,L,xn之间的点x处的函数值时，就需用一个较简单的且满足条件

(xi)=yi，i=0,1,2,L,n

的函数 (x)来代替f(x)，这一方法被称为插值法，其中 (x)被称为插值函数，f(x)被称为被插值函数，点xi被称为插值点．

当这一较简单的函数 (x)被选择为多项式时，则称为代数插值．当这一较简单的函数 (x)被选择为三角函数时，则称为三角插值．下面主要来介绍代数插值．

(1) 线性插值（也称一次插值）

已知函数表

x y

x0 y0

x1 y1

去构造一个线性多项式 1(x)，使其满足 1xi=yi，i=0,1． 这样的一次插值多项式具有如形式

1(x)=y0l0(x)+y1l1(x) 其中 l0(x)=

x x1x x0

，l1(x)=，被称为线性插值的基函数．误差为R(x)=f(x) 1(x)=Oh2，

x0 x1x1 x0

()

其中h=x1 x0为步长．

(2) 二次插值

已知函数表

x y

x0 y0

x1 y1

x2 y2

去构造一个二次多项式 2(x)，使其满足 2xi=yi，i=0,1,2．这样的二次插值多项式为

2(x)=y0l0(x)+y1l1(x)+y2l2(x) 其中 l0(x)=

(x x1)(x x2)，lx=(x x0)(x x2)，lx=(x x0)(x x1)，被称为二次

()x xx x()x xx x

x0 x1x0 x21101222021插值的基函数．误差为R(x)=f(x) 2(x)=O(h3)，其中h为步长．

(3) n次插值

已知函数表

x y

x0 y0

x1 y1

x2 y2

…… ……

xn yn

去构造一个n次多项式 nx)，使其满足 nxi)=yi，i=0,1,2,L,n．

n次插值多项式为

n(x)=其中 li(x)=

∑

i=0

n

yili(x)

(x x0)(x x1)L(x xi 1)(x xi+1)L(x xn)，i=0,1,2,L,n，被称为n次插值的基xi x0xi x1Lxi xi 1xi xi+1Lxi xn函数．误差为R(x)=f(x) n(x)=O(hn+1)，其中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

和n次插值函数的函数文件 (\chap1\interpolyp.m)

function y=interpolyp(x1,y1,x,n) a=interpolyb(x1,x,n); y=a*y1';

说明当n=1和n=2时，分别可以得到一次插值函数和二次插值函数．

例 在某实验中要测某短时光束的发光强度随时间的变化，已知在实验中测得的一组数据如下，t为测数据的时间，y为反映光束强度的指标，求t=2.25秒处的光的强度指标． 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\figure1_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];

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')

该程序的运行结果不仅可以显示出t=2.25秒处的光强度为 0.7781，同时可以画出光强度随时间的演

第二节 拟 合（Fitting）

正如上一节所讲，插值主要用于获取一个给定的离散集合的局部逼近．如果要获取一个给定的离散集合的整体行为，我们则采用拟合．也就是说，给定一组离散数据

x

x1 y1

x2 y2

n

x3

…… ……

xn

y

y3 yn

要求一个函数y=p(x，按照最小二乘原理，使得

2

[()]px y ∑iii=1

达到最小值，这一方法被称为拟合法，其中p(x)被称为拟合曲线．

值得注意的是拟合曲线y=p(x)的确定可根据实验或实际情况来选择．通常我们考虑多项式拟合，

即选取拟合曲线p(x)为m次多项式p(x)≡pm(x)=使得函数

∑

k=0

m

akxk，则问题变为求ak(k=0,1,2,L,m)，

f(a0,a1,Lam)=∑

i=1

n

m ∑akxik yi k=0

2

达到最小值．确定ak(k=0,1,2,L,m)的方法是，让

f

=0, aj

于是得到

j=0,1,2,L,m

∑

k=0

m

n

nk+j

∑xi ak=∑yixij,

i=0 i=1

j=0,1,2,L,m

这是一个m+1元线性方程组．

特别地，当m=1时，即线性拟合，此时拟合直线p1(x)=a1x+a0由系数a1和a0完全确定，它们满足如下方程组

nn

n2

a1∑xi+a0∑xi=∑xiyi i=1i=1i=1

nn

a

1∑xi+a0n=∑yi i=1 i=1

通过引入向量矩阵

x11 y1 x1 y A= 2 和 b= 2

MMM x1 y n n

我们可以将该方程改写为如下的矩阵形式

a1 T

AA a =Ab

0

T

由此我们就可以得到拟合直线的系数为

a1 T =AA a 0

()

1

ATb．

下面是为求拟合直线的系数a1和a0所编写的MatLab函数文件 (\chap1\fitlinec.m)

function yy=fitlinec(x,y) 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

x=x(1:n);y=y(1:n);

x=reshape(x,n,1);y=reshape(y,n,1); A=[x ones(n,1)];

B=A'*A; b=A'*b; yy=B\b; yy=yy';

和函数文件 (\chap1\fitlinep.m)

function p1=fitlinep(x,y,t) a=fitlinec(x,y); p1=a*[t;1];

一般地，我们引入矩阵

x1mx1m 1Lx11 y1 m m 1

xnLx21 x y2

和 =bA= 2

MMOMM M

y xmxm 1Lx1

n nn n

则拟合曲线pm(x)的系数am,am 1,L,a1,a0所满足的方程组可以改写为

am a

ATA m 1 =ATb

M a 0

由此我们就可以得到拟合曲线的系数为

am 1La0)=ATAATb．

下面是为求拟合曲线pm(x)的系数am,am 1,L,a1,a0所编写的MatLab函数文件 (\chap1\fitpolyc.m)．

T

(am

()

1

function 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<=m

warning('The number of the DATA should be greater than m'); return; end

x=x(1:n);y=y(1:n);

x=reshape(x,n,1);y=reshape(y,n,1); A=ones(n,1); for i=1:m

A=[x.^i A];

b=y; B=A'*A; b=A'*b; yy=B\b; yy=yy';

和拟合曲线pm(x)的函数文件(\chap1\fitpolyp.m) function pm=fitpolyp(x,y,t,m) a=fitpolyc(x,y,m); T=1;

for i=1:m T=[t^i T]; end

pm=a*T';

特别地，当m=1时，可以归结为直线拟合的情况．

例 某实验中测得一组数据，其值如下

x = 1 2 3 4 5 y =1.3 1.8 2.2 2.9 3.6．

那么m次多项式拟合曲线与数据点的关系可以用如下程序的计算结果及图象来表述 (\chap1\ figure1_2_1.m)

x=[1 2 3 4 5];

y=[1.3 1.8 2.2 2.9 3.6]; m=2; n=10;

h=(max(x)-min(x))/n; for i=1:n+1

t(i)=min(x)+(i-1)*h;

pm(i)=fitpolyp(x,y,t(i),m); end

plot(t,pm) hold on

plot(x,y,'b*')

运行结果给出拟合曲线与数据点的关系，可以用如下图象来表述．

图1.2.1 拟合曲线与数据点的关系

第三节 数值微分 (Numerical Differentiation)

据

若函数f(x)由表格给出，求f(x)在节点上的微商，称这样的问题为数值微分，即给定一组离散数

x y

x0 y0

x1 y1

x2 y2

…… ……

xn yn

求函数f(x)在节点xi(i=0,1,2,L,n)上的微商（其中f(xi)=yi）．

采用插值多项式来代替f(x)，去求函数在节点xi(i=0,1,2,L,n)处的微商．有三种思路：(1) 先插值后计算微商；(2) 利用Taylor展开；(3) 插值与微商同时进行（如样条函数）．

下面我们利用Taylor展开来给出数值微分中的对称三点式、前向两点式、后向两点式和对称五点

式的公式以及二阶导数的差分格式．

已知

x y

… …

x 2h f 2

x h f 1

x f0

x+h f1

x+2hf2

… …

求函数f(x)在节点x处的微商（近似值）．

由f(x)在x点处的Taylor展开

2 x)(+

3

x)(+

f(x+ x)=f(x)+ x f'(x)

知

2!

f′′(x)

3!

f′′′(x)+LL

h2h3h4(4)

f±1=f(x±h)=f0±h f'(x)+f′′(x)±f′′′(x)+f(x)+Οh5， (1.3.1)

2624

4h32h4(4)2

f±2=f(x±2h)=f0±2h f'(x)+2hf′′(x)±f′′′(x)+f(x)+Οh5， (1.3.2)

33

由此可以得到下面几个数值微分公式：

(1) 对称三点式

由(2.1)得到

h3

f1 f 1=2h f'(x)+f′′′(x)+Οh5，

3

即

f1 f 1h2

f′(x)=f′′′(x)+Οh4， 2h6

取

f f 1

f′(x)≈1，其截断误差为Οh2． (1.3.3)

2h

其中公式(2.3)被称为对称三点式公式．另外需要说明的是，若f(x)在[x h,x+h]上是二次多项式，

()

()

()

()

()

则该估算式是精确的．实际上，该估算式是x h,x,x+h三点插值多项式的微商值．

(2) 前向两点式

由(2.1)中f1的表达式

h2h3h4

f1 f0=h f'(x)+f′′(x)+f′′′(x)+f

2624

知

(4)

(x)+Ο(h5)，

f1 f0=h f'(x)+Οh2，

由此得到

()

f1 f0

，其截断误差为Ο(h)． (1.3.4) h

这里公式(2.4)被称为前向两点式公式．若f(x)在[x,x+h]上是一次多项式，则该估算式是精确的．实

f′(x)≈

际上，该估算式是x,x+h两点一次插值多项式的微商值．需要强调的是公式(2.4)的精确度比对称三点式低一个量级．

(3) 后向两点式

由(2.1)中f 1的表达式

f 1

知

h2h3h4

f0= h f'(x)+f′′(x) f′′′(x)+f

2624

f 1 f0= h f'(x)+Οh2，

(4)

(x)+Ο(h5)，

()

由此得到

f′(x)≈

f0 f 1

，其中截断误差为Ο(h)． (1.3.5) h

这里公式(2.5)被称为后向两点式公式．若f(x)在[x h,x]上是一次多项式，则该估算式是精确的．实际上，该估算式是x h,x两点一次插值多项式的微商值．同样地，公式(2.5)的精确度比对称三点式低一个量级，与前向两点式是同一量级．

(4) 对称五点式

为了得到更高的精确度，同时利用(2.1)和(2.2)，则有

h3

f1 f 1=2h f'(x)+f′′′(x)+Οh5，

38h3

f2 f 2=4h f'(x)+f′′′(x)+Οh5，

3

从中消去f′′′(x)，得到

()

()

8(f1 f 1) (f2 f 2)=12h f'(x)+Οh5，

由此得到

()

f′(x)≈

其中截断误差为Οh4，我们称公式(2.6)被称为对称五点式公式．若f(x)在[x 2h,x+2h]上是四次多项式，则该估算式是精确的．实际上，该估算式是x 2h,x h,x,x+h,x+2h五点四次插值的微商值．注意该公式的精确度比对称三点式高两个量级．

(5) 二阶微商的数值微分

()

1

f 2 8f 1+8f1 f2)， (1.3.6) (12h

利用(2.1)知

f1+f 1

由此得到

h4

=2f0+h f′′(x)+f

12

2

(4)

(x)+Ο(h5)，

f′′(x)≈

(6) 五点二阶微商的数值微分

利用(2.1)和(2.2)知

f1+f 1 2f0

h2

，其中截断误差为Οh2． (1.3.7)

()

从中消去f

(4)

(x)，我们得到

h4(4)

f1+f 1=2f0+h f′′(x)+f(x)+Οh6，

12

4h4(4)2

f2+f 2=2f0+4h f′′(x)+f(x)+Οh6，

3

2

()

()

f′′(x)≈

其中截断误差为Οh

()．

4

1

[16(f1+f 1) (f2+f 2) 30f0]， (1.3.8) 12h2

需要说明的是以上各公式均被称为差分公式．

第四节 数值积分 (Numerical Integration)

本节目的是介绍计算积分

∫f(x)dx的近似方法，主要介绍梯形公式和Simpson公式．

a

b

1.梯形公式

利用定积分的几何意义可以导出如下梯形公式

bhf(x)dx=f(a)+2f(a+h)+2f(a+2h)+L+2f(a+(n 1)h)+f(b)+Oh2， a2

∫

[

]()

其中h=

b a

． n

2.抛物线（Simpson）公式

计算积分

∫f(x)dx的Simpson公式为

a

b

∫

其中h=

b

a

h

f(x)dx=[f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+2f(a+4h)

3

+L+4f(a+2(n 1)h)+2f(a+(2n 1)h)+f(b)]+Oh3，

()

公式的推导：首先将区间[a,b]分成n等份，每个小区间为[xi,xi+1]，i=0,1,2,L,n 1，这样有

b a

． 2n

∫

对于每个小区间[xi,xi+1]上的积分

b

a

f(x)dx=∑

i=0

n 1

∫

xi+1

xi

f(x)dx．

∫

xi+1

xi

f(x)dx，我们用二次插值函数来代替其中的被积函数，于是得

到

∫

b

a

f(x)dx=∑

i=0

n 1

∫

xi+1

xi

f(x)dx

=

1

(xi+1 xi)f(xi)+4fxi+xi+1+f(xi+1) 6h

=[f(a)+4f(a+h)+2f(a+2h)+4f(a+3h)+2f(a+4h) 3

+L+4f(a+2(n 1)h)+2f(a+(2n 1)h)+f(b)]，

())

其中积分

∫

xi+1

xi

f(x)dx可以利用Maple7来实现，其中程序为：

> 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])*((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

(4fi + 1/2 + fi + 1 + fi)(xi xi + 1) 6

i+1

这里fi=f(xi),fi+1/2=f

xi+xi+1

),f

=f(xi+1)．

下面我们举几个利用MatLab计算数值积分的例子．

例1：计算积分

∫

1

exdx．

显然，该积分的精确值为

∫

1

exdx=e 1≈1.718281828．下面我们分别利用梯形公式和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)*h); end

int=h/2*(int+integrand1(b)); format long int 输出结果为

int =1.71828182989095，

其中被积函数的函数文件程序为：

function y=integrand1(x) y=exp(x);

利用Simpson公式的MatLab程序：(\chap1\simpson1.m)

% Using Simpson rule to calculte the integral. a=0; b=1; n=100;

h=(b-a)/(2*n);

int=integrand1(a); for i=1:n i;

int=int+4*integrand1(a+(2*i-1)*h); end

for j=1:n-1

int=int+2*integrand1(a+2*j*h); end

int=h/3*(int+integrand1(b)); int 输出结果为

int =1.71828182846501．

从上述的程序和计算结果可以看到，利用梯形公式循环10000次所得的计算结果，如果采用Simpson公式的话，只需循环100次即可．这说明Simpson公式具有比梯形公式更好地精度，当然这也提高了计算的速度．

例 2：计算积分

∫

1

x 2/3(1 x)

1/3

dx．

这是一个瑕积分的例子，有两个奇点x=0和x=1，一般的处理方法是将积分区域分为若干个部分，然后对每个积分作不同的变量替换，消除掉奇点．就本例来说，具体的处理形式为

∫ =∫

1

18

x 2/3(1 x)x

2/3

1/3

dx dx+7x

12

12

(1 x)

1/3

1

2/3

(1 x)

1/3

dx+1x 2/3(1 x)

78

78

1/3

dx

=

∫(

31 x3

)

1/3

dx+∫3x1 x3

()

2/3

dx+1x 2/3(1 x)

8

1/3

dx．

作为练习．

例3：对于脉宽大于5ps的光脉冲，描述其在单模光纤内传输的方程是非线性薛定谔方程

β2 2Aα A2

i= iA+ γ|A|A， (1.4.1) 2 z22 T

方程(1.4.1)右端的三项分其中A是脉冲包络的慢变振幅，T=t z/vg是脉冲随vg移动的时间度量．

别对应于光脉冲在光纤中传输时的吸收效应、色散效应和非线性效应．根据入射脉冲的初始宽度T0和峰值功率P0，决定脉冲在光纤内演变过程中是色散还是非线性效应起主要作用．

我们首先将方程(1.4.1)无量纲化，为此引入归一化振幅U和归一化初始脉宽T0

A(z,τ)=P0e

αz/2

Tt z/vg

U(z,τ)，τ==， (1.4.2)

T0T0

式中P0是入射脉冲的峰值功率，指数因子代表光纤的损耗．这样方程(1.4.1)就变为

Usgn(β2) 2Ue αz2i= |U|U， (1.4.3) z2LD τ2LNL

式中sgn(β2)=±1，它根据群速度色散参量β2的符号来确定，而

1

(1.4.4)

|β2|γP0

分别被称为色散长度LD和非线性长度LNL，它们分别表示沿光纤长L方向脉冲演化过程的长度量．根据色散长度LD、非线性长度LNL和光纤长度L的相对大小，脉冲演变可分为下面讨论的四种不同的传

LD=

，LNL=

输区域．

(1) 在L<<LD，L<<LNL的情况

当光纤长度L<<LD，L<<LNL时，色散和非线性效应都不起重要作用，此时方程(1.4.3)右端的两项可被忽略．于是我们可以得到，U(z,τ)=U(0,τ)，即脉冲在传输过程中保持其形状．在这个区域，光纤不起太重要的作用，只起传输光脉冲的作用，除了由于吸收引起的脉冲能量的降低，因而此区域对光通信系统是有益的．

(2) 在L≈LD，L<<LNL的情况

当光纤长度L≈LD，L<<LNL时，方程(1.4.3)右端的最后一项可被忽略，此时在脉冲的演变过中群速度色散起主要作用，而非线性效应相对较弱．由方程(1.4.4)可以看出，当光纤和脉冲参量满足下述关系时，适用于以色散为主的区域

T02

γP0T02LD=<<1 (1.4.5) |β2|LNL

(3) 在L<<LD，L≈LNL的情况

当光纤长度L<<LD，L≈LNL时，方程(1.4.3)右端的第一项较非线性项可以忽略．在这种情况

下，光纤中脉冲的演变过程自相位调制起主要作用，它将导致脉冲的频谱展宽．由方程(1.4.4)可以看出，当光纤和脉冲参量满足

γP0T02LD=>>1 (1.4.6) |β2|LNL

时，适用于以非线性为主的区域．

(4) 在L≥LD，L≥LNL的情况

当光纤长度L≥LD，L≥LNL时，脉冲在光纤内传输过程中，色散和非线性效应将共同起作用．群速度色散和自相位调制效应的互作用与群速度色散或自相位调制效应单独起作用相比较有不同的表现．在反常色散区（β2<0），光纤能维持光孤子．在正常色散区（β2>0），群速度色散和自相位调制可用来进行脉冲压缩．当群速度色散和自相位调制效应都重要时，方程(1.4.3)对理解光纤中脉冲的演变是很有帮助的．

在本例中，我们将考虑在L≈LD，L<<LNL的情况，此时方程(1.4.3)或方程(1.4.1)右端的最后一项可被忽略，可以写成

Aβ2 2A

， (1.4.7) i=2

2 T z

方程(1.4.7)只包括了最低阶的群速度色散β2，尽管大多数实际情况下这一项的影响是主要的，但有时也需要考虑更高阶的色散，如三阶色散β3．此时方程(1.4.7)可以改写为

β3 3A Aβ2 2A

(1.4.8) i=+i

2 T26 T3 z

此方程能利用傅立叶方法求解，其解的形式为

β331+∞~ β22

A(z,T)=Aωiωz+iωz iωT(0,)exp dω (1.4.9)

262π∫ ∞

~

其中入射场A(0,T)的傅立叶变换A(0,ω)为

+∞~

A(0,ω)=∫A(0,T)exp(iωT)dT． (1.4.10)

∞

若入射场A(0,T)是确定的，方程(1.4.9)就能用来研究高阶色散效应．

下面我们来计算初始高斯脉冲的演化，其中A(0,T)具有如下形式

A(0,T)=exp( T2/2T02)． (1.4.11)

三阶色散β3可以被忽略，此时可以通过直接积分求得(1.4.9)在大多数情况下群速度色散β2是主要的，

的表达式．然而，当β2=0时，高阶色散效应β3起主要作用，此时我们很难求得积分(1.4.9)的表达式．这样为了看到高阶色散的作用，数值计算积分(1.4.9)就显得非常重要．为此我们引入与高阶色散

′ 有关的三阶色散长度LD

′=T03/|β3|． (1.4.12) LD

这样我们就可以将积分(1.4.10)和(1.4.9)写成如下形式

~

A(0,ω)=T0∫

=T0

和

+∞

∞+∞

∫

∞

1 T 2 T T

exp expiωTd 0 2 T0 T0 T0

1

exp τ2 exp(iωT0τ)dτ≡T0(0,ωT0)，

2

+∞

A(z,T)=

β3T 3z dωT0 (0,ωT)expi(ωT)iωT0003∫ ∞ 6T0 T0 13z1+∞T = i (0, )exp i ∫ d ， (1.4.13) ∞′2π6LTD0 1

2π

其中

(0, )=∫

+∞

∞

1

exp τ2 exp(i τ)dτ． (1.4.14)

2

并假定β3>0．下面是我们编写的程序，其中计算积分(1.4.14)的函数M文件为(\chap1\gsint.m) function q=gsint(omega) b=10;n=200;ht=2*b/n; t=-b:ht:b; int1=0; for j=1:n

int1=int1+ht/2*(GS(t(j))*exp(i*omega*t(j))+GS(t(j+1))*exp(i*omega*t(j+1))); end

q=int1;

而计算积分(1.4.13)的函数M文件为(\chap1\gsint1.m)

function q=gsint1(Z,t) b=10;n=200;hw=2*b/n; w=-b:hw:b; int=0; for j=1:n

int=int+hw/2*(gsint(w(j))*exp(i/6*w(j)^3*Z-i*w(j)*t)... +gsint(w(j+1))*exp(i/6*w(j+1)^3*Z-i*w(j+1)*t)); end

q=int*conj(int)/(2*pi)^2;

最后可以得到运行程序(\chap1\gsplot.m)

clear all;

a=5;n=5;hz=a/n;

b=10,m=200,ht=2*b/m; z=0:hz:a; t=-b:ht:b; for j=1:n+1 j

for k=1:m+1

q(j,k)=gsint1(z(j),t(k)); end end

mesh(q)

这里初始入射高斯脉冲的函数M文件为(\chap1\GS.m)

function q=GS(t) q=exp(-t^2/2);

数值结果如图1.4.1和图1.4.2所示．程序中我们采用了梯形公式．另外需要说明的是，由于这里的积分是无穷积分，一般的处理方法是将积分区间限制在一个有限的、大的范围内进行数值积分．

从图1.4.1和图1.4.2可以看出，三阶色散会引起脉冲形状的畸变，会在其中一个脉冲沿附近形成非对称的振荡结构．正像图1.4.1和图1.4.2所示，当β3>0时，振荡出现在脉冲的后沿．事实上，当β3<0时，振荡会出现在脉冲的前沿．

练习题：用Maple来推导梯形公式． 练习题：计算积分

∫

+∞

∞

e xdx．程序见\chap1\trapezoid2.m, simpson2.m 和integrand2.m．

2

第五节 求 根

在许多实际问题中常常会遇到求解非线性方程和非线性方程组的问题．本节目的是介绍求解代数方程f(x)=0的根的近似方法．

求解方程的根，通常会遇到有两种情况，一种是要求给定范围内的某个根，而根的粗略位置可以从问题的物理背景或实际应用得知．另一种方法是求出方程的全部根，而根的数目和位置事先不知道，这在解超越方程时是比较困难的．

1.区间对分法

设函数f(x)在区间[a,b]上连续，且在(a,b)内只有一个零点．为方便起见，这里不妨假设

f(a) f(b)<0，且f(a)<0,f(b)>0．根据连续函数的性质，方程f(x)=0在(a,b)区间内必有一根，

区间(a,b)被称为方程的有根区间．

区间对分法是求解方程f(x)=0在(a,b)内的根有效方法，具体思路如下：首先取区间[a,b]之中点

a+ba+b a+b a+b a+b

，并计算f 即为方程的解．否则，若f ，如果恰好f =0，则 <0，

2 2 2 22

a+b

则令此中点为a1，且将b改写为b1；若f >0，则令此中点为b1，而将a改写为a1．这样我们

2

得到一个新的有根区间(a1,b1)，其长度是原来有根区间(a,b)的一半．如此下去，重复上述过程n次，我们可以得到一个有根区间序列

(a,b),(a1,b1),(a2,b2),L,(an,bn)

其中后一个有根区间的长度是前一个有根区间长度的一半．这样最后一个有根区间(an,bn)的长度为

(a,b)的1/2n，此时可取(an,bn)的中点作为方程f(x)=0在(a,b)区间内的根一个近似值，即

a+bn

xn≈n，

2

它和精确解的误差小于

b a

． n+1

2

区间对分法的优点是简单易实现，收敛速度与比值为求实的单根，不能求重根．

1

的等比级数相同，它的局限性是只能用于2

3

例1：利用区间对分法求解方程x 2x 5=0在区间(2,3)内的根． 下面是利用区间对分法求解该方程的程序(\chap1\ex1_5_1.m) clear all

a=2; b=3; x=(a+b)/2; Tol=10^(-15); j=1

while (b-a)/2>Tol

j=j+1

if f1(x)==0 return

else if f1(x)<0 a=x; x=(a+b)/2; else

b=x; x=(a+b)/2; end end

end

format long x

运行结果为2.09455148154233．

2.迭代法

对于给定的方程f(x)=0，可以将其改写为等价形式x= (x)，例如可令 (x)=x+f(x)．由该等价形式，我们可以非常方便地写出如下形式的迭代格式

xn+1= (xn),n=1,2,3,L．

n→∞

这样的话，给定一个初始猜测x1，就可以得到一串数列{xn}．如果数列xn收敛，即limxn=ξ，那么有

ξ=limxn+1=lim (xn)= (ξ)

n→∞

n→∞

这说明ξ是方程f(x)=0的根，这样，我们就可以将xn视为f(x)=0的近似解．总之，若函数f(x)连续且数列xn收敛，则由迭代格式得到的数列xn即为f(x)=0的近似解．

迭代过程的几何意义：把方程f(x)=0求解的问题，改写为求数列{xn}的极限ξ，实际上是把问题转化为求直线y=x和曲线y= (x)的交点，其横坐标就是ξ，如下图所示．从图中可以看出，用同样的办法从x1出发构造出的数列{xn}，有的数列离ξ越来越近，它被称为收敛的数列，如图1.5.1；有的离ξ越来越远，它被称为是发散的数列，如图1.5.2．那么，在什么条件下，迭代数列是收敛的？可以证明：在直线y=x和曲线y= (x)的交点邻近，如果函数 (x)的Lipschitz常数小于1，那么迭代数列是收敛．通常函数 (x)的Lipschitz常数小于1的条件可以用条件′(x<1来代替．

下面我们考虑一个例子：

图1.5.1 迭代法收敛示意图

图1.5.2 迭代法发散示意图

例2：用迭代法求解方程x+2x+10x 20=0． 我们将该方程写成两种等价形式

32

x=x+x3+2x2+10x 20 和 x=20/x2+2x+10．

32

分析发现，由第一种形式给出的 (x)=x+x+2x+10x 20在与直线y=x的交点邻近不满足迭代

()

数列收敛的条件，而由第二种形式得到的 (x)=20/x+2x+10在与直线y=x的交点邻近满足条件 ′(x)<1，所以我们采用第二种形式．由此得到如下迭代格式

2

xn+1=20/xn+2xn+10

(

2

)

()

并取初始值为x1=1．下面是由MatLab编写的程序（\chap1\ex1_5_2.m）

clear all

x1=1;

x=f2(x1);

Tol=10^(-15); i=1;

while abs(x-x1)>Tol i=i+1 x1=x;

x=f2(x1); end x

运行结果为1.36880810782137．从这个例子我们可以看到，迭代方法的收敛性与函数 (x)的选取有关，所以当采用迭代方法求根时，选取一个满足条件 ′(x)<1的函数 (x)是至关重要的．

3.搜索法

值得注意的是，以上所介绍的方法适用于求解方程在指定位置的一个根，当我们要同时求解多个根时，一般地来说我们应该采用搜索法．利用搜索法求根时的基本要求是已知根的大体位置（通常利用分析的方法或其它实验手段获取）．下面我们介绍搜索法的基本思路．

首先在根的存在范围内任取一x0，不妨设x0的值比根小（当x0的值比根大时也可同样考虑）；然后以一正步长h增加x0的值为x1=x0+h，并比较函数值f(x0)和f(x1)．若f(x0) f(x1)<0，则将

步长h减半倒退回来，即x1

h

x1；否则在x0+h的基础上再增加步长h，重复上述过程，如图所2

示．随着循环次数的增加，这种方法产生的值必收敛到根．值得注意的是，搜索法具有方向性，搜索方向的选择是至关重要的，如果搜索方向选择错误，则会出现求不出解的情况．

图1.5.3 搜索法示意图

2

例3：利用搜索法数值计算方程x 5=0的根．

容易得到，该方程的精确解为x=±5=2.236067978． 下面是利用搜索法求解该方程的程序(\chap1\searching.m) clear all; t=10^(-6);

h=0.5; x=0;

while abs(h)>t fo=f(x); x=x+h;

if fo*f(x)>0 else

x=x-h; h=h/2; end end x

其中函数文件为(\chap1\f.m)

function p=f(x) p=x^2-5;

如果要求全部的两个根，则以上程序可以改写为(\chap1\fzeros.m)

% This is a program which evaluate the zeros of

% the function f(x) by using the searching aigorithm

clear all;

t=10^(-15); a=-5;

for i=1:2 h=0.1; a=a+h; fo=f(a);

while abs(h)>t a=a+h;

if fo*f(a)>0 else

a=a-h; h=h/2; end end

x(i)=a; end x'

运行结果为 -2.23606797749979

2.23606797749979

例4：对于脉宽T0<1ps的超短脉冲，需要考虑包括高阶色散和高阶非线性效应在内的高阶非线性薛定谔方程

本文来源：https://www.bwwdw.com/article/xblm.html