偏微分方程数值解

这一组装过程通过命令行eassempde来实现。

约束处理，求解方程组

对于Neumann条件，由于是自然边界条件，??上不需要满足任何约束条件，因此可立即得到线性代数方程组：KU=F 一旦形成K和F，在MATLAB环境下立即可解出节点近似解的向量u。如果是Dirichlet边界条件，还需要对边界节点进行约束处理，然后再解方程组。 4．2 抛物型方程

下面说明抛物型方程如何简化成椭圆型方程来求解。这是通过PDE Toolbox函数parabolic完成的。

?u????c?u??au?f,in?,考虑?t （4-7）

d初值为

u(x,0)?u0(x),x??. （4-8）

边界条件类似椭圆边值问题的提法，这里仅讨论Neumann条件：

n?(c?u)?qu?g,??. （4-9）

对于热传导方程（4-7）可以写成

?C

?u???(k?u)?h(u?u?)?f,?t （4-10）

表示热量向环境扩散，其中是热源。

?是密度，C是比热，k是导热系数，h是薄层传热系数，u?是环境温度，f

如果系数与时间无关，方程就是标准的椭圆型方程： -??(c?u)?au?f.

对?作三角形网格剖分，对于任意给定t?0,PDE的解按有限元法的基底可以展开成

u(x,t)??ui(t)?i(x).i （4-11）

将展开式带入方程（4-7），两边乘以试验函数可得

?i，

然后在?上积分，利用Green公式和边界条件（4-9），

??d?j?idxi??dui(t)??(???j?(c??i)?a?j?idx??q?j?ids)vi(t)dxdti??? （4-12）

??f?jdx??g?jds,?j.??M上式可写成大型的线性稀疏的常微分方程组： 这就是所谓的线性半离散化方法。

再解ODE（4-13）的初值，初值为

dU?KU?F.dt （4-13）

Ui(0)?u0(xi).

得每一个节点

xi、任一时刻t的ODE解。这里K和F是远边界条件下椭圆型方程

???(c?u)?au?f,in?

的刚度矩阵和荷载向量。M是质量矩阵。

当边界条件是与时间有关的Dirichlet条件hu=r时，F项则包括h和r的时间倒数，它可以用有限差分求解。常微分方程组是病态的，这时需作显式时间积分。

由于稳定性要求时间间隔很小，而隐式解由于每一时间段都要求解椭圆型方程，从而求解非常缓慢。常微分方程组的数值积分，可以由MATLAB中suite函数完成。对于这类问题它是有效的。 4．3 双曲型方程

类似于抛物型方程的有限元法，可以解双曲型问题：

?2ud2???(c?u)?au?f,in?. ?t u(x,0)?u0(x),?u(x,0)?v0(x),x??,?t边界条件同上。

初值为

?2u?c?u?02?t特别地，对于波动方程，波速为c。

对区域?作三角形网格剖分，与抛物型方程处理方法一样，可以得到二阶常微分方程组：

d2VM2?KV?F,dt

Vi(0)?u0(xi),dVi(0)?V0(xi),?i.dt

初值为

K和M是刚度矩阵和质量矩阵。

PDE Toolbox中提供的求解双曲型方程的命令函数是hyperbolic。 4．4 特征值方程

在PDE Tooltox中求解的基本特征值问题是 ????c?u??au??du，

其中?是未知复数。在固体力学中方程描述薄膜震动的问题，在量子力学中?描述势阱a(x)中有界定态的能级。数值解包括方程的离散和代数特征值问题的求解。首先考虑离散化。按有限元基底将u展开，两边同乘基函数，再在?上作积分，可以得到广义特征值方程：

KU=?MU，

其中对应于右端项的质量矩阵的元素为

Mi,j??d(x)?j(x)?i(x)dx.?

利用命令函数assema（单元积分贡献的装配）生成。在通常情况下，当函数d（x）为正时，质量矩阵M为正定对称矩阵。同样，当c（x）是正的而且在Dirichlet边界条件下，刚度矩阵K也是正定的。

广义特征值问题：

KU=?MU.

利用Arnoldi算法进行移位和求逆矩阵,直到所有的特征值都落在用户事先确定的区间.在此,求解的详细过程略去.

PDE Toolbox中提供的求解特征值问题的命令函数是pdeeig。