Matlab微分方程高效解法：一维泊松方程求解

在狄利克莱边界条件下，考虑一维泊松问题：

u′′=f(x),-1<x<1,u(±1)=0 （5-23）

其中，f(x)为已知函数，u(x)为待求函数。将横轴上的区间[-1,1]离散化为向量x=(x₀,x₁,…,x_N)^T，相应地，u(x)被离散化为向量u=(u₀,u1,…,uN)T。则式（5^-23）可以写为矩阵形式：

D_N²u=f(x),u₀=u_N=0（5-24）

式（5-24）代表了含有N+1个方程的方程组，其中未知数有N-1个：u₁,u₂,…,u_N-1，比方程个数少2个。实际上，只需要N-1个方程就能求解这些未知数，不妨忽略边界上的方程，即：删去矩阵D_N²的首尾行以及向量f(x)的首尾元素。又因为u₀=u_N=0，矩阵D_N²的首尾列与其相乘的结果也是0，所以删去矩阵D_N²的首尾列和向量u的首尾元素。于是就得到了在狄利克莱边界条件u(±1)=0下的2阶切比雪夫求导矩阵——将矩阵D_N²的首尾行、首尾列删除后的N-1阶方阵，如图5-6所示，此矩阵将在后续内容中反复用到。

图5-6 在狄利克莱边界条件u(±1)=0下修改切比雪夫求导矩阵

如果这里用“~”表示删除矩阵首尾行和首尾列、删除向量首尾元素的操作，那么求解式（5-23），只需：

注意，必须先对D_N平方，再删除其首尾行和首尾列，这个顺序不能颠倒。此外，通过式（5-25）得到N-1维向量后，一定要在其首尾补0。

通过在切比雪夫求导矩阵上乘以缩放因子，上述方法可以从区间[-1,1]推广到任意区间，以下面的泊松方程为例：

u′′=x²-x,-2＜x＜2,u(±2)=0 （5-26）它的精确解为u(x)=x⁴/12-x³/6+2x/3-4/3。计算数值解并与精确解比较的代码如下：

程序5-5

程序输出结果如图5-7所示。N=10时，精确解与数值解吻合得很好，最大误差在10^-15数量级，而且边界处的误差最小。

图5-7 左图：精确解（曲线）与数值解（点）的最大误差在10^-15数量级。右图：误差分布

如果泊松方程的狄利克莱边界条件更具一般性，比如：(www.xing528.com)

那么，有两种方法解决这一问题。

方法1：

式（5-27）的解可写为u(x)=u₀(x)+c₁x+c₂的形式，其中u₀(x)是满足式（5-26）的特殊解。选取合适的c₁、c₂使多项式c₁x+c₂满足式（5-27）的边界条件，即：(c₁x+c₂)|_x=2=2，(c₁x+c₂)|_x=-2=-1，得到c₁=0.75，c₂=0.5。毫无疑问，u₀(x)+0.75x+0.5必是式（5-27）的解。

方法2：

如图5-8所示，将矩阵D_N²的首尾行分别修改为100…0和0…001，以使得向量u和向量f(x)的首尾元素相等。然后将向量f(x)的首尾元素分别改为u₀和u_N（函数u(x)在边界的取值），最后通过下式计算向量u：

图5-8 在更具一般性的狄利克莱边界条件下修改切比雪夫求导矩阵

其中，矩阵和向量上的“^”代表对矩阵、向量进行上述修改。在方法2中，不但利用了N-1个方程求解未知数u₁,u₂,…,u_N-1，还额外加了两个方程（对应于矩阵D_N²的首尾行），以确保u的首尾元素必满足边界条件。

用方法1和方法2求解式（5-27）并与精确解u(x)=x⁴/12-x³/6+17x/12-5/6比较的代码如下：

程序5-6

程序输出结果如图5-9所示，在N=10的情况下，方法1的计算结果和方法2的计算结果均与解析解吻合较好，误差分别在10^-15和10^-16数量级。此外，方法1在各处的误差要普遍高于方法2在各处的误差。

图5-9 左图：方法1的计算结果（○）、方法2的计算结果（+）、解析解（曲线）；右图：方法1的误差（○）、方法2的误差（+）