模拟物理-09 第七章 椭圆型偏微分方程

x
• 运行这个程序我们可以看到，尽管对解初始猜测 与真实解相去很远，但是迭代仍然收敛。求得的 解与解析结果有良好的一致。
• 收敛的速度依赖于松弛参数ω。一般选择它大于1。
最佳值可以由经验方式决定。
• 上述松弛方法可以用于到二维（甚至三维）问题。
• 二维的算法是，迭代过程中
• 把这个算法相继应用于网格上的每一点，比方说 按从上到下的顺序扫描各行，而在同一行中从左 到有扫描，就收敛到所需的解。
椭圆型偏微分方程
偏微分方程
• 任何一种物理现象，若其物理量随空间变化或者 同时随着空间和时间变化，那么对这种现象的描 述中一定包含有偏微分方程。
• 这些现象时多种多样的，比如扩散、电磁场的分
布、电磁波、流体力学和量子力学（Schrödinger
波）。
• 除了最简单的情形之外，这些方程是不能用解析 方法求解的，必须用数值方法才能得到定量的结 果。 • 在典型的数值处理方法中，我们选取自变量（即 空间和时间）的许多离散的点，计算其上的因变 量（例如电位，温度等）。通过适当的离散化， 偏微分方程就转化为一大组差分方程。
• 把已知量（边界）移动到等号右边
• 写成矩阵
• 对于任何一种实际情况，方程的数量是巨大的。 • 比如当N=50则方程个数约为2500。 • 因此用直接求M的逆矩阵的方法来求解线性方程组是 不实际的。 • 由于微分方程的离散化近似只包含相邻的点，矩阵的 大部分元素为零，矩阵是稀疏的。 • 求解这种方程，存在几种效率很高的迭代方法。
双曲型方程
• 双曲型方程含有符号相反的二阶微商。例如描述 一根绷紧的弦的波动方程。
• 本章讨论适用于椭圆型方程的一些数值方法。
• 第七章讨论抛物型方程。 • 双曲型方程常常可以用相似的方法处理。
• 在本章我们将讨论适用于椭圆型方程的一
些数值方法。
• 为了讨论具体起见，我们将考虑关于二维 空间(x,y)内的场ϕ的椭圆型方程边值问题
写出方程的矩阵形式，然后用矩阵运算的方法求
解。
• 这里我们要讨论一个非常有效的特殊方法。
• 首先我们把差分方程改写 • 这个方程并不是显式的方程，因为我们并不知道方程 右边的各个ϕ的值。
• 我们的基本想法（Gauss-Seidel迭代）是：先猜测一个
初始解，然后对格子进行系统的扫描（比方说从左到
右），相继地把每一点上的ϕ换成一个由上述方程给
• 在我们开始讨论迭代方法时，我们先考虑一个类 似的但是更简单的一维边值问题，然后再回到二 维情形上来。
• 同我们讨论的二维问题相似的一维边值问题可以
写为
• 边条件为给定ϕ(0)和ϕ(1)的值。
• 差分方程为
• 针对这个离散方程，我们前面讨论过使用向前积
分或向后积分的方法。
• 我们也可以把所有方程中的变量组成一个列向量，
• 大部分物理上重要的偏微分方程是二阶的，可以 分为三种类型：抛物型、椭圆型和双曲型。
抛物型方程
• 抛物型方程对一个变量只有一阶微商，但
是含有对其它变量的二阶微商。
• 其例子是扩散方程和含时间的Schrödinger 方程。它们对时间的微商是一阶的，对空
间的微商是二阶的。
椭圆型方程
• 椭圆型方程包含对每个自变量的二阶微商， 并且当方程中所有的项集中到方程的一边 时，每项的符号都相同。 • 这类方程包括静电场电位的Poisson方程， 拉普拉斯方程，赫姆霍兹方程，和不含时 间的Schrödinger方程。它们都含有两个或更 多个空间变量。
• 这些差分方程在原则上可以用前一章讨论的直接 矩阵方法求解，但是涉及的矩阵非常大。矩阵的 维数与格点的数量相当，常常大于几千。庞大的 矩阵使得这种方法是不实际的。 • 幸好，原始的微分方程的局域性使得生成的差分 方程是稀疏的，即有关的矩阵的大部分元素为零。 对于这样的矩阵，使用各种迭代方法来求逆和使 矩阵对角化可以有很高的效率。这些方法是本章 的主题。
• 在(i, N − 1)点上离散化的方程，可以改写为
• 偏微分方程的上述离散化近似等价于关于区域的 内点上的未知量ϕ值的一个线性方程组。 • 用矩阵记号，它可以写为
• ϕ是用所有内点排列成的列向量。 M是出现在方
程组中的矩阵，S包含了离散化方程右边的两项。
例子
• 对于一个单位正方形，每条边分成4等份。格点编 号是0,1,2,3,4 • 内点是：(1,1), (1,2), (1,3), (2,1), (2,2), (2,3), (3,1), (3,2), (3,3) • 把内点代入方程：
• 下面的程序实现了松弛算法。其中使用的初始猜
测解是ϕi = 0.
program main implicit none integer, parameter :: N=20 real :: phi(0:N),S(0:N),x,phip real :: omega=1.0 real :: h=1.0/N integer :: i,it
离散化
• 我们的第一步是要把方程改写为一个适合进行数 值处理的形式。 • 为此，我们定义一个网格，覆盖(x, y)平面内我们 感兴趣的区域。 • 为了方便起见，我们取格子间隔h为均匀的，并且 在两个方向上相等，使得单位正方形被(N + 1) × (N + 1)个格点覆盖。
• 这些格点可以用指标(i, j)编号，每个指标从
! 初始化
do i=0,N
x=h*i
s(i)=h*h*12*x*x
phi(i)=0 ! 初始猜测值 enddo
! 迭代循环 do it=1,500 do i=1,N-1 phip=(phi(i-1)+phi(i+1)+S(i))/2 phi(i)=(1-omega)*phi(i)+omega*phip enddo enddo
Hale Waihona Puke Baidu
0到N，使得第(i, j)点的坐标为(xi = ih, yj = jh)。
• 如果我们再定义ϕ ij = ϕ(xi , yi )，并类似地定 义Sij ，那么对每个方向上的二阶微商应用
三点差分近似，就可以把方程(6.1)近似为
• 现在我们必须讨论边界条件在何处进入离散化的方程。 • 如果边界和网格的格点重合，可以直接进入后面的讨 论。 • 如果边界不规则，无法与格点重合，格点只能粗略地
• “超松弛”对应于ω > 1，而“低松弛”则意味着
ω < 1。ω的取值范围是0到2。
• 作为这个松弛方法的一个例子，我们考虑下述边 值问题： • 方程为
• 边条件为ϕ(0) = ϕ(1) = 0。 • 它的精确解是������ ������ = ������(1 − ������ 3 ).
• 首先我们写出这个方程的差分公式
• 方程的求解还需要边界条件。 • 我们将取Dirichlet型边界条件，（稳定场方 程的第一类边界条件）。 • 即在(x, y) 平面内的某一根很大的闭合曲线 上规定了ϕ的值。 • （也许还在其中的某些附加曲线上规定了ϕ 的值）。
• 为了方便，可以把边界闭合曲线取为单位
正方形。
• 别的边值问题，可以用非常相似的办法处 理。
描述边界的形状。那么可以使用非均匀的网格，在边
界处增加更多的格点，减小格点与边界的差别。
• 不论怎样，边界条件提供了在格点的一个子集上的ϕij
的值。
• 在远离一条边界的一点上，边界条件不直接进入 方程。 • 考虑和一条边界相邻的一点，例如(i, N − 1)点上的
离散化的方程，可以改写为
• 如果边界条件是第二类边界条件（Neumann边界 条件），比方说∂ ϕ /∂y = g(x)， • 那么边界条件可以近似为下述离散形式
• 虽然这不是椭圆型方程的最普遍的形式，但是 它已经代表了物理中许多的情况下的方程形式。 • 例如，在一个静电学问题中， ϕ是位势，S与 电荷密度有关； • 在一个定态热扩散问题中， ϕ是温度，S则是 局部的热量产生或损失的速率。 • 这里的考虑可以直接推广到其他情形，比如三 维空间的情形，或者扩散常数或介电常数随空 间变化的情形。
出的值。
• 注意在使用方程的时候，右边要用ϕi ± 1 的
最新的值。
• 多次重复这一扫描过程，就可以把对ϕ 的一 个初始猜测“松弛”到正确的解。
• 对于这个迭代过程，我们可以作一个推广，使得 松弛过程的每一步ϕi换成它的老值和方程给出的 值的一个线性组合
• 其中ω是一个参数，可以调节它来控制松弛速度：
! 输出 open(20,file=’data.txt’)
do i=0,N
x=i*h
write(20,*)x,phi(i),x*(1-x**3)
enddo close(20) endprogram
0.5 0.4 0.3

0.2 0.1 0.0 0.0
松弛方法 解析解
0.2
0.4
0.6
0.8
1.0