计算物理课件

物理问题



连续体(气体、液体、固体)和场(电磁场、温度场) 状态用关于时间和空间的函数描述：u = u(x, y, z, t) 系统遵循特定的物理规律，即 u 的变化满足特定的(偏微 分)方程 稳定过程：泊松方程 静电场：电场的散度正比电流密度，等于势场u的梯度

数学方程的建立


静磁场：类似于静电场
√
一维波动方程(3/3)

例：求解一维波动方程
2u 2 2u = 0, 0 x l , 0 < t < tmax 2 -v 2 x t u ( x, t ) = x(1 - x) u ( x, t ) t =0 = sin x, t t =0 u (0, t ) = u (l , t ) = 0
√
二维扩散方程(2/3)

交替方向隐式差分格式(ADI 格式)



由 k 时的 u 值求 k1 时的 u 值，由 k1 时的 u 值求 k2 时的 u 值 uy"在 k1/2 时用 k 的中心差商(反向的)，在 k3/2 时用 k2 的中心差商(正向的)——交替方向的 要求解线性方程组——隐式的 误差 O(h2) O(t 2)，无条件稳定的 线性方程组形式上与一维扩散方程的相同，解法也相同
√
一维扩散方程(2/3)

隐式六点差分格式(C-N 格式)


由 k 时的 u 值求 k1 时的 u 值 要求解线性方程组——隐式的 涉及六个结点的 u 值——六点 误差 O(h2) O(t 2)，无条件稳定的

边界条件的差分格式
√
一维扩散方程(3/3)

差分方程组及其求解
√
二维扩散方程(1/3)

差分方程组的特点



方程个数等于内点数，每条方程最多含5个未知项 系数矩阵是稀疏和带状的 跌代法求解：同步法、异步法和逐次超松弛法

同步法


用第 k 步的 u 值，代入方程的右边，计算得到新的 u 值， 用在第 k1 步，…；直到新旧的 u 差值小于设定的误差 特点：需要两套内存(存放 u )，收敛慢

二维扩散方程
2u 2u u , 0 < x < l x , 0 < y < l y , 0 < t < t max D( 2 2 ) = y t x u ( x, y,0) = u0 ( x, y ) u a1u b1 n = c1 ( y, t ), x = 0 u a u b = c2 ( y, t ), x = l x 2 2 n a3u b3 u = c3 ( x, t ), y = 0 n u a4u b4 = c4 ( x, t ), y = l y n


√
迭代解法(5/6)

例：求解拉普拉斯方程
2u 2u 2 2 = 0, 0 < x < 4, 0 < y < 3 x y u x =0 = y ( y - 3), u x =5 = 0 x , u y =4 = 0 u y =0 = sin 4
计算物理
有限差分方法
http://125.217.162.13/lesson/ComputationalPhysics
有限差分方法

物理问题和数学方程 有限差分原理 矩形区域中的泊松方程 迭代解法 非矩形区域中的泊松方程 一维扩散方程 二维扩散方程 一维波动方程
√
物理问题和数学方程(1/5)
√
一维波动方程(2/3)

显式差分格式


由 k-1 和 k 时的 u 值，直接求 k1 时的 u 值，不必解方程 组——显式的 误差 O(h2) O(t 2)，当 P<1 时收敛和稳定

初值条件的差分格式

边界条件的差分格式
2hc1 4b1u2 - b1u3 u = 1 2ha1 3b1 2hc2 4b2u N -1 - b2u N - 2 u N = 2ha2 3b2


抛物型(B2-4AC=0)：一维扩散方程 2u u = , { x, t , A = 1, B = C = 0} 2 x t 双曲型(B2-4AC>0)：一维波动方程 2u 2u - 2 = 0, { x, t , A = 1, B = 0, C = -1} 2 x t
√
物理问题和数学方程(2/5)

输运过程：扩散方程 扩散：流体由于不均匀而发生扩散，扩散密度正比于 密度的梯度，同时满足质量守恒

热传导：类似于扩散

振动传播：波动方程 交变电磁场

机械振动
√
物理问题和数学方程(3/5)

方程的分类
二阶线性偏微分方程的一般形式 2u 2u 2u u u A 2 B C 2 = F ( , , u, , ) 2 椭圆型(B -4AC<0)：二维泊松方程 2u 2u 2 = f ( x, y), { x, y, A = C = 1, B = 0} 2 x y
√
物理问题和数学方程(4/5)

边界条件和初始条件
物理：过程的具体特征与初始状态和边界(受外界影响)有关 数学：偏微分方程有无限个解，定解需要初始和边界条件 边界条件 第一类 u = u0 (rb , t ), 是区域D的边界，u0是已知函数， rb是上的位矢 第二类 u = q0 (rb , t ), q0是已知函数， n是的外法线 n u 当 u 是电磁场的势，则 代表场强 n u 当 u 是密度/温度/位移，则 是流量/热流/应力 n 当 q0=0 时称为第二类齐次边界条件

一维波动方程
2u 2 2u = f ( x, t ), 0 x l , 0 < t < t max 2 -c 2 x t u ( x, t ) = R ( x), u ( x, t ) = R ( x) 1 2 t =0 t t =0 a u b u = c , x = 0 1 1 1 n u = c2 , x = l a2u b2 n
√
有限差分原理(3/3)

误差为O(h2)的差商公式 dui - ui 2 4ui 1 - 3ui 一阶向前： = O(h 2 ) dx 2h dui ui - 2 - 4ui -1 3ui 一阶向后： = O(h 2 ) dx 2h dui ui 1 - ui -1 一阶中心： = O(h 2 ) dx 2h d 2ui ui 1 - 2ui ui -1 2 二阶中心： 2 = O ( h ) 2 dx h 收敛性：当步长 h0 时，差分方程的解是否收敛于微分 方程的解 稳定性：误差 Du 在运算过程中是否失控，即累计误差是 否无限增加


求解方法

u = f ( r ) 初始瞬间 u 在各处对时间的变化率： 2 t t =0
大多数没有解析解，只能数值求解 求解方法：变分解法、有限元法、有限差分法，边界元法
√
有限差分原理(1/3)

差商格式

基础：用差商代替微商(离散化)
√
有限差分原理(2/3)

误差为O(h)的差商公式

异步法

计算第 k1 步的 ui,j 时， ui-1,j 和 ui,j-1 已经知道 特点：需要一套内存(存放 u )，收敛较快
√
迭代解法(2/6)

逐次超松弛法


将第 k 步的 u 值，与异步法中计算得到的 u 值，加权平均， 用在第 k1 步 特点：加权因子 w 对跌代次数的影响很大。经验公式
√
uC 26.786
非矩形区域中的泊松方程(1/3)

二维圆形区域(极坐标系(r,j))

泊松方程的差分格式

圆心(i =1)的差分格式：以上格式不适合， 利用直角坐标系的五点差分格式

周期性条件
Baidu Nhomakorabea
边界条件
√
非矩形区域中的泊松方程(2/3)

例：同轴线两导线(内导线半径 a，外导线半径 b)间为均 匀介质(相对介电常数 e)，电位差为 V0，求电势分布 2 2 u 1 u 1 u 2 u= 2 2 = 0, a r b 2 r r r r j ua = 0, ub = V0
√
二维扩散方程(3/3)

边界条件的差分格式

虚点：i=0, i=N1, j=0, j=M1

回代交替方向隐式差分格式，得到边界的差分格式 当 t=1,3,5,..., 时，依次选择 x，沿 y 方向计算方程组 当 t=2,4,6,..., 时，依次选择 y，沿 x 方向计算

√
一维波动方程(1/3)
√
作业

求解泊松方程。取 h=0.2, w=1.25, e=10-4。并与解析 解 u=xey 比较
2 2 u u 2u = 2 2 = xe y , 0 < x < 2, 0 < y < 1 x y u x =0 = 0, u x =2 = 2e y , u y =0 = x, u y =1 = xe
√
迭代解法(6/6)

例：求解拉普拉斯方程，边界条 件如右图所示


取 h=5，有三个内点：uA, uB, uC 转化为线性方程组


方程的解 uA = 1.786, uB = 7.143 , uC = 26.786 收敛性
h uA uB 5 1.786 7.143 5/2 1.289 6.019 26.289 5/4 1.144 5.632 26.144 5/8 1.107 5.525 26.107 5/16 1.097 5.498 26.097 5/32 1.095 5.491 26.095

差分格式
√
非矩形区域中的泊松方程(3/3)

轴对称区域(柱坐标系(r,f,z)(r,z))

泊松方程的差分格式

对称轴(i =1)的差分格式：以上格式不适合

其它形状区域/边界

网格边界/近似插值实际边界 有限元方法/边界元方法
√
一维扩散方程(1/3)

一维扩散方程
2u ( x, t ) u ( x, t ) = , 0 x l , 0 t t max D 2 x t u ( x, t ) t = 0 = u 0 ( x ) u a u b = c1 , x = 0 1 1 n a2u b2 u = c2 , x = l n

求解一维扩散方程。取 a1=b1=a2=-b2=1, c1=c2=0, l=1, tmax=10, D=0.1, h=0.1, t=10-4。并与解析解 u=e x0.1t 比较
√
迭代解法(3/6)

矩形区域的第二和三类边界条件

当 a 和 b 是 x, y 的函数时，应 a=a(xi, yj) 和 b=b(xi, yj) 对第二边界条件，令 a=0
√
迭代解法(4/6)

不规则区域

第一类边界(不对称网格方法)

第二类边界 结点在边界上 结点不在边界上：过结点 P 向边界作垂线， 交于 P' 点，以 P 代替 P' 第三类边界 前两类边界条件的组合
√
物理问题和数学方程(5/5)

第三类

u (a0u b0 ) = c0 (rb , t ), a0、b0和c0是已知函数 n 例：热传导，系统通过表面与外界交换热量：表面 u 热流 正比于表面温度 u 与外界温度 u0 之差，即 n u = k (u - u0 ) n 初始条件 初始瞬间待求函数 u 在各处的值： u = f ( r ) 1 t =0
√

差商格式的收敛性和稳定性


矩形区域中的泊松方程(1/1)

五点差商格式

二维泊松方程 收敛性和稳定性：在给定边界条件下具 有唯一解，当 h0 时趋于解析解 方程( f(x,y)=0 ) 2 2 u u 2 u= 2 2 =0 x y 五点差商格式


拉普拉斯方程


√
迭代解法(1/6)