有限差分法

有限差分法是以差分原理为基础的一种数值计算法。它用各离散点上函数的差商来近似替代该点的偏导数，把要解的边值问题转化为一组相应的差分方程。然后，解出差分方程组（线性代数方程组）在各离散点上的函数值，便得边值问题的数值解。

现以二维等步长差分格式为例，说明有限差分法的原理和方法步骤。

1.区域离散化，作网格剖分

图1⁃4⁃1 二维等步长正方形网络

如图1⁃4⁃1所示，用平行于坐标轴的两组直线族将地下划分成正方形网格，相邻两坐标线的距离为h，则任一点的x、z坐标为

x=ih（i=0，1，2，…，M）

z=kh（k=0，1，2，…，N）

每个正方形为一单元，其边长h称为步长，网格的交点称为节点。任一节点的坐标（x，z）可表示为（ih，kh），或简化为（i，k），用阶梯状折线代替原来的曲线段。在边界线以内的节点称为内节点，边界上的节点称为边界节点。

2.微分方程离散化，构组差分方程

某一内节点（i，k）处的电位为U（i，k），由于h很小，可将节点（i，k）四周的电位在节点处展成泰勒级数：

地电场与电法勘探

地电场与电法勘探

地电场与电法勘探

地电场与电法勘探

式中Ux，Uxx，……和Uz，Uzz，……分别表示U对x和z的一阶导数、二阶导数等。将前两个式子相加，并且忽略h的四次项与更高次项，经整理可得：

地电场与电法勘探

同理得：

地电场与电法勘探

将上述Uxx和Uzz代入含源分区均匀岩石中位函数U所满足的微分方程（1⁃4⁃16）的第二式，即得二维函数U（x，z）的差分方程：

U（i+1，k）+U（i，k-1）+U（i-1，k）+U（i，k+1）-4U（i，k）=h2f（1⁃4⁃18）

对于无源分区均匀介质，位函数 U（x，z）所满足的微分方程（1⁃4⁃17）的差分方程为

U（i+1，k）+U（i，k-1）+U（i-1，k）+U（i，k+1）-4U（i，k）=0（1⁃4⁃19）

3.线性方程组的形成与求解

对于边界节点，其相应的差分方程可根据边界条件给出。全部结点所建立差分方程（1⁃4⁃18）和（1⁃4⁃19）的总和可分别写成以下矩阵形式：

〔A〕·｛U｝=｛F｝（1⁃4⁃20）

和

〔A〕·｛U｝=0（1⁃4⁃21）

〔A〕是方程组的系数矩阵，它是与电阻率分布有关的函数；｛U｝是电位U的列向量，其分量为所有节点上的电位；｛F｝是常向量。当给定电阻率分布及边界条件后，解线性方程（1⁃4⁃20）和（1⁃4⁃21），便可求得电位的空间分布。

电位｛U｝值的计算精度与步长h的大小有很大关系。一般说来，网格划分越细，即h值越小，｛U｝值与理论值就越接近。但是此时节点数目也急剧增加，因而所需的计算机内存和计算时间也就会增大。解决计算速度与精度这一矛盾的较好方法是采用变步长，即在近区将网格分得密些，远区影响较小可分得稀些。

如何使用matlab，用有限差分求解偏微分方程？

求解思路：把偏微分方程离散化，采用合适的差分方法，将复杂的方程简化成简单的线性方程组，最后求解线性方程组，得到其数值解。

现以一维扩散方程为例，说明其计算过程。

第一步，根据条件，建立边界条件和初始条件，即

g0=@(t)zeros(size(t))

g1=g0%边界条件

eta=@(x)sin(pi*x)%初始条件

第二步，设置网格数，即

n=101%网格数

m=101%网格数

第三步，设置步长，即

h=0.01%步长

k=0.01%步长

第四步，设置t和x的初始值，即

t0=0%t的初始值

x0=0%x的初始值

第五步，确定扩散系数，即

K=1/pi^2

第六步，自定义Crank-Nicolson差分格式解函数

[t,x,U]=diffusion_sol1(h,k,t0,x0,n,m,eta,g0,g1,K)

第七步，绘制偏微分方程解的曲面，即

surf(t,x,U)

最后，运行程序得到一维扩散方程数值解的曲面图

---