有限元分析（FEA）是个什么东东

J Pan

航空工程师

一、有限元能干什么
二、偏微分方程如何求解（有限元发明以前）
三、有限元法的基本思路
四、有限元法的数学基础——降维
五、如何获得“弱形式”的解
六、二维、三维有限元计算

一、有限元能干什么

通常能通过第二牛顿定律能得出这样的方程：

m\ddot x(t)+c \dot x(t) +kx(t)=f(t)

我们把这种方程称之为 常微分方程 ，即 自变量只和时间 t 有关系，和空间位置没关系 ，换句话说，常微分方程描述的是单质点的变化规律，如：某个物体在重力作用下做自由落体运动，下落距离随时间变化的规律；火箭在发动机推动下在空间飞行，飞行的轨道等等。常微分方程一般是把研究对象当成一个质点或者刚体，研究整体的运动规律。

值得高兴的是，线性常微分方程相对比较好解决，我们可以通过 傅里叶变换 （ J Pan：傅里叶变换后面的到底有什么小秘密）和 拉普拉斯变换 （ J Pan：从另一个角度看拉普拉斯变换）将微分方程变成代数方程，进而得出准确的 解析解 。

事实证明，对于繁复纷杂的大自然，只用常微分方程是不够的，因为有些研究对象不能简化成质点，一个典型的例子就是琴弦：

琴弦是一个柔性体，在弹拨的时候每个点振动都是不一样的，一根均匀的弦，其自由振动的数学方程是什么呢？假定 u(x,t) 表示 x 点在 t 时刻的位移，如果取琴弦中一个微元进行力学分析，就可以得到琴弦应遵守的方程为：

\frac{\partial^2u}{\partial t^2}-a^2 \frac{\partial^2u}{\partial x^2}=0

可以发现，这个方程的 自变量不仅仅是时间，还有空间坐标 ，我们把这种方程称之为偏微分方程。也就是说，偏微分方程能描述的连续体的各个点随时间连续变化的情况，本质上就是一种“场”的描述了。

严格来说， 自然界的各种现象，都需要用偏微分方程来描述 ，只不过有些情况下为简单起见，我们忽略了位置项带来的影响，如火箭发射的时候，在强烈的振动下壳体会发生变形的，理论上要完整描述火箭的状态，是需要用偏微分的，但是实际上我们还是把它简化成一个质点来描述，因为这样可以简化很多工作。

当然，这个琴弦除了满足上面的方程，还有其它约束，比如说琴弦的两端要固定在支架上，琴弦内部与外界是通过这个边界联系起来的，边界的变化是会影响琴弦的波动的，我们把这些称之为 边界条件 ，比对于琴弦，边界条件一般为：

u\big|_{x=0}=0 ； u\big|_{x=l}=0 ；

除了边界，琴弦的波动显然还和它的“历史”有关。两根同样的弦，一根在薄刀背的敲击下发生的声音就比较刺耳，另一根在手指的弹拨下就比较和谐，两根弦由于初始时刻的振动情况不一样，后来振动下发出来的声音就不一样，我们把初始时刻的状态称之为 初始条件 ： u\big|_{t=0}=\varphi(x) ； \frac{\partial u }{\partial t}\bigg|_{t=0}=\psi(x) ；

通常边界条件（含边界条件和初始条件）有三类，一类是变量本身的约束，如 u\big|_{t=0}=\varphi(x) ，称之为狄利克雷边界条件；一类是变量的导数有约束，如 \frac{\partial u }{\partial t}\bigg|_{t=0}=\psi(x) ，称之为诺依曼边界条件；还有第三类就是前两者的混合，成为罗宾边界条件。

遗憾的是，对于偏微分方程，除了少数极简单的情况下能计算出解析解（形式还很复杂），绝大多数情况下我们还无法很好的求解——有限元就是在这种情况下应运而生，简要来说， 有限元法就是为了求解偏微分方程！

二、偏微分方程如何求解（有限元发明以前）

前面我们说了，对于大多数偏微分方程而言，我们是没有办法求得解析解的，只能退而求其次，寻求个近似解，即便是这个近似解的过程也是相当不容易的，方法也很简单、粗暴—— 试凑法 ！比如我们现在我们假定一个位移函数，其中包含若干个待定系数：

u=u_0+\sum_{m}{A_mu_m}

当然这些函数也不是随意假定的，它需要满足位移边界条件，这个靠 u_0 来满足， u_m 主要用来模拟非边界条件的点，它们在边界处取值应该为零，否则和 u_0 叠加后会破坏位移边界条件。 A_m 为待定系数。

系统在外力的作用下，发生位移，产生变形。位移可以是各种各样的，但必须满足位移的边界条件。但是，满足位移边界条件的位移位移有无穷多组，我们应该选哪一组呢？——那就是当 总势能取极小值 时的位移，这种现象，我们称之为 最小作用量原理 ，这是自然界中一个普遍成立的原理。

实际我们是怎么做的呢？我们知道，不同的 A_m 组合就会获得不同的位移，这样我们就可以列出总势能的表达式，其中只包含待定系数 A_m ，由前面我们说的最小势能原理， A_m 只有一组是真实的解，那就是当总势能变分 \Pi 为零的时候，怎么做到的呢？很简单，将 \Pi 分别对 A_m 求导取零就行了，这个时候偏微分方程组就变成线性方程组了，这种解偏微分的方法称之为利兹法。

如果我们所取的位移不仅满足位移边界条件，而且根据它们求得的应力还满足应力边界条件(不要求满足平衡方程），这种方法称之为伽辽金方法，这种方法对位移函数的要求较高，但计算量小一些。

上面说的都太抽象，我们举个具体的例子，注意下面的例子不用看懂，只是为了演示传统求解偏微分方程是多么费劲！

我们举一个利兹法的例子，可能大家会更容易理解一些：两端简支的等截面梁，受均匀分布载荷 q(x) 作用如图所示，试求解梁的挠度 u(x) 。

首先构造位移试函数（说是构造，其实就是猜），

u(x)=\sum_m\varphi_msin\frac{m\pi x}{l}

当然也不是瞎猜，这个试函数要满足位移边界条件

u\big|_{x=0}=0 ； u\big|_{x=l}=0 ；

则总的势能为（边看形式复杂，就是按公式）：

E_t=\frac{EJ}{2}\int_0^l (\frac{d^2 u(x)}{dx^2})^2dx-\int_0^l qu(x)dx

整理得到

E_t=\frac{EJ\pi^4}{4l^3}\sum_m m^4\varphi_m^2-\frac{2ql}{\pi}\sum_{m=1,3,5,\hdots}\frac{\varphi_m}{m}

根据 \frac{\partial E_t}{\partial \varphi_m}=0 ，得到 \varphi_m=\frac{4ql^4}{EJ\pi^5m^5} \quad m=1,3,5,\hdots

所以 u(x)=\frac{4ql^4}{EJ\pi^5}\sum_{m=1,3,5,\hdots}\frac{1}{m^5}sin\frac{m\pi x}{l}

再次强调，这个例子看思路即可！

三、有限元法的基本思路

根据第二节的描述，我们知道传统解决偏微方程的方法（如利兹法），采用的试函数 \varphi_i(x) 非常复杂，而且是在全域定义的，工程上的物体（零组件或者系统）及其边界条件都太复杂了，我们几乎不能够构建一个很好的基于整个弹性体（全域）的位移函数。

那怎么办呢？——还记得圆周率是怎么计算的吗？其中有一种几何方法是这样的，就是用正多边形等效圆形，就像切西瓜一样，将圆切成有限个等腰三角形，示意图如下：

每个等腰三角形的面积 S_i=\frac{1}{2}R^2sin\theta_i

由所有的等腰三角形组成的正多边形的面积为： S_N=\sum_1^N S_i=\frac{1}{2}R^2Nsin(\frac{2\pi}{N})

我们知道， \lim_{x \rightarrow 0}{\frac{sinx}{x}=1} ，所以当 N\rightarrow \infty 时， \lim_{N \rightarrow \infty}{Nsin(\frac{2\pi}{N}) }=2\pi 正多边形就趋向于圆形了，此时可获得圆的面积为 S=\pi R^2 。也就是说我们可以采用“ 化整为零 ”的思想，将复杂不易求解的东西切成简单的容易计算的几何形状来进行等效。

也就是说，当整体曲面或者曲线较为复杂时，我们就可以把它切开，切开的每一个规则的小块称之为单元，单元与单元之间通过节点联系起来，整个弹性体就被划分成了有限个单元，简称有限元。对一个规则的单元（子域）再假设位移函数就简单多了，比如可以采用线性的试函数 a+bx (x\in [x_i,x_{i+1}]) ，试函数的要求也不高（连续性阶次较低），缺点就是计算量变大了，不过我们现在有计算机了，这都不是事。

简而言之， 有限元分析就是“化质为量”，就是将解微分方程这个难题（质）转化为大量的数值计算（量）。

参考文献：

[1] 曾攀-有限元分析基础教程

四、有限元法的数学基础——降维

很多童鞋都读过刘慈欣的《三体》，里面提到了一个很火的概念，叫做“ 降维打击 ”，什么意思呢？就是将攻击目标本身所处的空间维度降低，致使目标无法在低维度的空间中生存从而毁灭目标。其实“降维”这种做法我们早就使用了，比如地图就是典型的降维的结果，因为我们关心的是经纬度、不关心每块地区的海拔高度，这样就可以将三维的地球信息用二维的地图表示。

有限元分析其实也就是降维分析，将“无限维”降低到“有限维”进行计算。以上都是抽象的描述，接下来我们看看在数学上怎么描述。举个简单的例子，在每个节点，我们选择线性分段函数作为试函数 \varphi_i(x) ，它的特点是在节点 x_i 处取值为1，在 x_j ( j \ne i )处取0，其他位置线性变化，形状曲线如下：

数学表达式如下：

\varphi_i(x)=\left\{ \begin{array}{lcl} (x-x_{i-1})/(x_i-x_{i-1}) && if \ x\in [x_{i-1},x_i]\\ (x_{i+1}-x)/(x_{i+1}-x_{i}) && if \ x\in [x_{i},x_{i+1}]\\ 0 && otherwise\\ \end{array} \right.

因为 \varphi_i(x) 长得像个帽子，所以一般称之为hat function。这些函数有个特点就是“作用范围”很小，它们只在节点周围的若干个单元上有值，其它地方一概为0。接下来要引入一个非常重要的概念： \varphi_i(x) 作为基函数，可以张成空间 V_h ， V_h 中任一函数都可以表示成：

v(x)= a_1\varphi_1(x)+a_2\varphi_2(x)+\hdots +a_n\varphi_n(x)=\sum_{i=0}^na_i\varphi_i(x)

什么意思，大家估计都玩过积木，基础模块只有几种，通过各种不同的空间组合，却能拼出各种房子、动物等等。

前面我们说的 \varphi_i(x) 就可以看成是积木中的基础模块，我们称之为基函数；积木搭建的各种房子、动物等，我们称之为基函数张成的空间 V_h ； \varphi_i(x) 的线性组合就是 V_h 中一个函数向量。这个点我们要多说一下，因为对于没有接触过泛函分析的人，确实不太容易理解。我们都很熟悉泰勒分析： f(x)=f(x_0)+\frac{f'(x_0)}{1!}(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+\hdots +\frac{f^n(x_0)}{n!}(x-x_0)^n+R_n(x)

这个方程啥意思呢？——多项式 1,x,x^2,x^3,\hdots, x^n,\hdots 为基函数组成一个空间，空间里的任意函数（足够光滑）都可以由这些基函数线性叠加而来。再比如傅里叶分解也是这个意思：任意周期函数（满足狄利克雷条件），都可由正余弦函数叠加而来，而系数就是被分解函数与基函数的内积，也就是被分解函数在基函数上的投影！

说到投影大家都很熟悉了，比如X光透射，就是光透过三维物体留下的二维影像！

那函数 f(x) 是不是也可以通过计算其在空间 V_h 的投影来计算投影后系数 a_i 呢？——当然是可以的，比如我们可以这么计算：

\int_I(f(x)-f_{prj}(x))v(x)dx=0, \quad \forall v(x)\in V_h

其中， f_{prj}(x) 为 f(x) 在空间 V_h 的投影， v(x) 为空间 V_h 任意函数向量。既然 f(x)-f_{prj}(x) 与 V_h 中所有的函数 v(x) 正交，则说明 f_{prj}(x) 是 f(x) 在空间 V_h 中的投影，具体请看下图：

由于 v(x) 的任意性，为方便起见，我们可以选择基函数 \varphi_i(x) ，于是可以得到：

\int_I(f(x)-f_{prj}(x))\varphi_i(x)dx=0

由于 f_{prj}(x) 在空间 V_h 中，我们可以假设

f_{prj}(x)=\sum_{j=0}^na_j\varphi_j(x)

所以

\begin{equation} \begin{aligned} \int_If_{prj}(x)\varphi_i(x)dx&=\int_I\left( \sum_{j=0}^na_j\varphi_j(x) \right)\varphi_i(x)dx\\ &=\sum_{j=0}^na_j\int_I\varphi_j(x)\varphi_i(x)dx \end{aligned} \end{equation}

这个方程看着太复杂，我们假设：

M_{ij}=\int_I\varphi_i(x)\varphi_j(x)dx, \quad i,j=0,1,\hdots,n

b_i=\int_If(x)\varphi_i(x)dx, \quad i,j=0,1,\hdots,n

于是方程

\int_I(f(x)-f_{prj}(x))\varphi_i(x)dx=0

就变成了

\sum_{j=0}^nM_{ij}a_j=b_i,\quad i=0,1,\hdots,n

注意，这是一个 (n+1)\times(n+1) 线性方程组， a_j 为 n+1 个未知量，我们通过解这个线性方程组，就可以获得基函数的系数 a_j 进而获得函数 f(x) 的近似。式中 b_i 称为负载向量，由于历史原因， M 一般称之为质量矩阵。

接下来我们就详细看一下这个线性方程组具体长什么样。我们先画出相邻的两个hat functions：

由前面的分析，采用辛普森积分法则，可以得到：

\begin{equation} \begin{aligned} M_{i,i} &=\int_I\varphi_i^2(x)dx\\ &=\int_{x_{i-1}}^{x_i}\varphi_i^2(x)dx+\int_{x_{i}}^{x_{i+1}}\varphi_i^2(x)dx\\ &=\frac{h_i}{3}+\frac{h_{i+1}}{3},\quad i=1,2,\hdots,n-1 \end{aligned} \end{equation}

其中 x_i-x_{i-1}=h_i ， x_{i+1}-x_i=h_{i+1} ， M_{0,0} 和 M_{n,n} 是half hat，因此需要单独计算 M_{0,0}=h_1/3 ， M_{n,n}=h_n/3 ,同样可以计算 M_{i+1,i} 如下：

\begin{equation} \begin{aligned} M_{i+1,i} &=\int_I\varphi_i(x)\varphi_{i+1}(x)dx\\ &=\int_{x_{i}}^{x_{i+1}}\varphi_i(x)\varphi_{i+1}(x)dx\\ &=\frac{h_{i+1}}{6} ,\quad i=0,1,2,\hdots,n \end{aligned} \end{equation}

同样可以得到： M_{i,i+1}=h_{i+1}/6 。

我们对质量矩阵进行组装，可以得到：

M=\begin{bmatrix} \frac{h_1}{3} & \frac{h_1}{6}\\ \frac{h_1}{6} & \frac{h_1}{3}+\frac{h_2}{3} &\frac{h_2}{6}\\ & \frac{h_2}{6} & \frac{h_2}{3}+\frac{h_3}{3} &\frac{h_3}{6}\\ &&\ddots & \ddots & \ddots\\ &&& \frac{h_{n-1}}{6} & \frac{h_{n-1}}{3}+\frac{h_n}{3} &\frac{h_n}{6}\\ &&&& \frac{h_n}{6} & \frac{h_n}{3} \end{bmatrix}\quad

MATLAB代码如下：

function M = MassMat1D(x)
n = length(x)-1; % number of subintervals
M = zeros(n+1,n+1); % allocate mass matrix
for i = 1:n % loop over subintervals
    h = x(i+1) - x(i); % interval length
    M(i,i) = M(i,i) + h/3; % add h/3 to M(i,i)
    M(i,i+1) = M(i,i+1) + h/6;
    M(i+1,i) = M(i+1,i) + h/6;
    M(i+1,i+1) = M(i+1,i+1) + h/3;
end

采用梯形积分法则，可以可以计算负载向量 b_i 如下：

\begin{equation} \begin{aligned} b_i &=\int_If(x)\varphi_i(x)dx\\ &=\int_{x_{i-1}}^{x_{i}}f(x)\varphi_i(x)dx+\int_{x_{i}}^{x_{i+1}}f(x)\varphi_i(x)dx\\ &=f(x_i)(h_i+h_{i+1})/2 \end{aligned} \end{equation}

完整矩阵形式为：

b=\begin{bmatrix} f(x_0)h_1/2\\ f(x_1)(h_1+h_2)/2\\ f(x_2)(h_2+h_3)/2\\ \vdots\\ f(x_{n-1}(h_{n-1}+h_n))/2\\ f(x_n)h_n/2 \end{bmatrix}\quad

MATLAB代码如下：

function b = LoadVec1D(x,f)
n = length(x)-1;
b = zeros(n+1,1);
for i = 1:n
    h = x(i+1) - x(i);
    b(i) = b(i) + f(x(i))*h/2;
    b(i+1) = b(i+1) + f(x(i+1))*h/2;
end

现在我们来验证一下算法是否正确，假定 f(x)=x\cdot sinx,x \in [0 ,1] ,代码如下：

function y = Foo(x)
y=x.*sin(x)

我们来计算一下投影函数和原函数分别长什么样，代码如下：

n = 5 ;% number of subintervals
h = 1/n; % mesh size
x = 0:h:1; % mesh
y = Foo(x);% calc f(x)
M = MassMat1D(x); % assemble mass
b = LoadVec1D(x,@Foo) % assemble load
Pf = M\b; % solve linear system
plot(x,Pf,x,y); % plot  projection

计算结果如下，可见投影函数基本和原函数一致，但是误差比较大，这是因为我们节点去太少了，减小 h 后重合度会更高。

参考文献：

[1] Larson, Mats G., and Fredrik Bengzon. "The finite element method: theory, implementation, and practice."

五、如何获得“弱形式”的解

前面我们计算了函数 f(x) 的基函数空间的投影情况，我们说有限元是为了解决偏微分方程求解的，那 f(x) 和偏微分方程又有什么关系呢？这就牵扯到微分方程的“弱形式”了！ 获得微分方程的“弱形式”解，可以说是有限元法最重要的概念之一 ，我们看一下是怎么一回事！假设我们微分方程具有下列形式：

-\frac{\partial ^2u}{\partial x^2}=f(x), \quad x\in I=(0,1)

边界条件为：

\frac{\partial u}{\partial x}\big|_{x=0}=\kappa_0(u|_{x=0}-g_0)

-\frac{\partial u}{\partial x}\big|_{x=1}=\kappa_1(u|_{x=1}-g_1)

其中 k(x)>0 和 f(x) 为给定的函数， \kappa_0\geq0 ， \kappa_1\geq0 ， g_0 和 g_1 是给定常数。具有给定热源的一维热传导方程就具备这种形式。微分方程描述的是微观状态，但是我们一般更容易掌握的是宏观的现象，所以呢，方程我们一般喜欢把微分方程转化成积分方程，比如对上式两边同时积分：

\int_0^1-\frac{\partial ^2u}{\partial x^2}dx=\int_0^1f(x)dx

乍一看，皆大欢喜，实际上有问题，问题在哪呢？——以上方程表示在整个计算域 [0 \quad 1] 中，方程的积分值相等，也就是 平均值相等 ，与原来要求处处相等相比，显得“太弱了”！那怎么办呢？——我们可以找一个试函数 v(x) （任意的），将原方程改写成如下方式：

\int_0^1-\frac{\partial ^2u}{\partial x^2}\cdot v(x)dx=\int_0^1f(x)\cdot v(x)dx

当然， v(x) 不能表现太糟糕，至少要保证两边积分都是收敛的。上述方程对于所有满足条件的 v(x) 都成立，这个要求太强了，显然当 -\frac{\partial ^2u}{\partial x^2}=f(x) 时这个方程是成立的，那是不是只有这种情况下才成立呢？——实际上就是的，这可以通过变分法证明。

下一个问题就来了，既然 v(x) 有千千万万，那我们到底选择哪一个？总不能都选择吧！——前面我们说了，有限元就是为了解决偏微分方程的求解问题，思路就是将复杂的空间切成很多小区域，每个区域上假设一个试函数test function（或者叫形函数shape function） \varphi_i(x) ，然后通过求取试函数的系数 a_i 来获得函数 f(x) 的近似。试函数有一个特点，它只在一个很小的区间上有数值，其他地方都为零。如果我们在每个点上都用试函数测试一遍，这样至少能保证所单元上都成立的！——这种将试函数和形函数选择为同一个函数的方法，是一个俄国人最先使用的，我们称之为伽辽金（Galerkin）法。

这样我们就得到了如下的方程：

\int_0^1-\frac{\partial ^2u_h}{\partial x^2}\cdot \varphi_i(x)dx=\int_0^1f(x)\cdot \varphi_i(x)dx

这就是所谓的原微分方程的 弱形式 ！之所以“弱”，是因为这个方程的解只能保证在每个单元上方程左右两边的平均值相等，而不是原微分方程要求的每一点都相等。注意此时的 u_h(x) 已经与原微分方程的解 u(x) 有了差异，因为我们现在是用有限维空间去近似了原无限维空间。

接下来我们要继续变形了，要稍微用点数学知识——不知道大家还记不得分部积分法：假设有函数 w 和 v ，则很容易得到 wv 乘积的微分为：

d(w\cdot v) =w dv+vdw+dw\cdot dv

dw \cdot dv 为二阶无穷小，可略去，所以可以得到：

wdv=d(wv)-vdw

两边积分得到： \int wdv=\int d(wv)-\int vdw=wv-\int vdw

现在我们令 w=\varphi_i(x) ， v=\frac{\partial u}{\partial x} ，则有

\begin{equation} \begin{aligned} \int_0^1-\frac{\partial ^2u_h}{\partial x^2}\cdot \varphi_i(x)dx &=-\int_0^1 \varphi_i(x)d\left( \frac{\partial u_h}{\partial x} \right)\\ &= \int_0^1 \frac{\partial u_h}{\partial x} \frac{\partial \varphi_i}{\partial x}dx-\left( \frac{\partial u_h}{\partial x} \cdot {\varphi_i(x)}\right)\bigg|_{x=0}^{{x=1}}\\ \end{aligned} \end{equation}

所以我们可以得到：

\int_0^1f(x)\cdot \varphi_i(x)dx=\int_0^1 \frac{\partial u_h}{\partial x} \frac{\partial \varphi_i}{\partial x}dx-\left( \frac{\partial u_h}{\partial x} \cdot {\varphi_i(x)}\right)\bigg|_{x=0}^{{x=1}}

通过这种变化，我们对 u(x) 的要求由2阶降到了1阶，将形函数 \varphi_i(x) 从0阶增加到了1阶，由于形函数是我们自己设计的，这就降低了对原微分方程解的光滑性要求。将边界条件代入上式，可得展开的微分方程弱形式：

\begin{equation} \begin{aligned} &\int_0^1 \frac{\partial u_h}{\partial x} \frac{\partial \varphi_i}{\partial x}dx +\kappa_1u_h(x)\varphi_i(x)\bigg|_{x=1} +\kappa_0u_h(x)\varphi_i(x)\bigg|_{x=0}\\ &=\int_0^1f(x)\cdot \varphi_i(x)dx+\kappa_1g_1\varphi_i(x)\bigg|_{x=1} +\kappa_0g_0\varphi_i(x)\bigg|_{x=0} \end{aligned} \end{equation}

我们可以假设

u_h(x)=\sum_{j=1}^{n-1}a_j\varphi_j(x)

代入展开的微分方程弱形式就可以得到：

(A+R)a=b+r

其中 A 和 R 为 (n+1) \times (n+1) 矩阵， b 和 r 为 n+1 向量，具体表达式如下：

A_{i,j}=\int_0^1\frac{\partial \varphi_j}{\partial x} \frac{\partial \varphi_i}{\partial x}dx

R_{i,j}=\kappa_1\varphi_i(x)\varphi_j(x)|_{x=1}+\kappa_0\varphi_i(x)\varphi_j(x)|_{x=0}

b_i=\int_0^1f(x)\varphi_i(x)dx

r_i=\kappa_1g_1\varphi_i(x)|_{x=1}+\kappa_0g_0\varphi_i(x)|_{x=0}

这样我们就将微分方程转化成了线性矩阵方程：

\sum_{j=1}^{n-1}(A_{i,j}+R_{i,j})a_j=b_i+r_i

其中 A 为 (n-1) \times (n-1) 矩阵，称之为刚度矩阵， b_i 为负载矩阵，和前面的定义一样。和前面的计算相似，可以获得 A 的完整矩阵形式如下：

A=\begin{bmatrix} \frac{1}{h_1} & \frac{-1}{h_1}\\ \frac{-1}{h_1} & \frac{1}{h_1}+\frac{1}{h_2} &\frac{-1}{h_2}\\ & \frac{-1}{h_1} & \frac{1}{h_2}+\frac{1}{h_3} &\frac{-1}{h_3}\\ &&\ddots & \ddots & \ddots\\ &&& \frac{-1}{h_{n-1}} & \frac{1}{h_{n-1}}+\frac{1}{h_n} &\frac{-1}{h_n} \\ &&&& \frac{-1}{h_n} & \frac{1}{h_n} \end{bmatrix}\quad

R_{i,j} 为边界条件带来的等效矩阵，显然只有 i=j=0 或者 i=j=n 时 R_{i,j} 才有取值，很容易求得 R_{0,0}=\kappa_0 ， R_{n,n}=\kappa_1 ，即

R=\begin{bmatrix} \kappa_0\\ &&&\\ &&&\kappa_1 \end{bmatrix}\quad

所以可以得到：

A+R=\begin{bmatrix} \frac{1}{h_1} & \frac{-1}{h_1}\\ \frac{-1}{h_1} & \frac{1}{h_1}+\frac{1}{h_2} &\frac{-1}{h_2}\\ & \frac{-1}{h_1} & \frac{1}{h_2}+\frac{1}{h_3} &\frac{-1}{h_3}\\ &&\ddots & \ddots & \ddots\\ &&& \frac{-1}{h_{n-1}} & \frac{1}{h_{n-1}}+\frac{1}{h_n} &\frac{-1}{h_n} \\ &&&& \frac{-1}{h_n} & \frac{1}{h_n} \end{bmatrix}+\begin{bmatrix} \kappa_0\\ &&&\\ &&&\kappa_1 \end{bmatrix}

MATLAB代码如下：

function A = StiffMat1D(x,kappa)
n = length(x)-1;
A = zeros(n+1,n+1);
for i = 1:n
    h = x(i+1) - x(i);
    A(i,i) = A(i,i) + 1/h; 
    A(i,i+1) = A(i,i+1) - 1/h;
    A(i+1,i) = A(i+1,i) - 1/h;
    A(i+1,i+1) = A(i+1,i+1) + 1/h;
A(1,1) = A(1,1) + kappa(1);
A(n+1,n+1) = A(n+1,n+1) + kappa(2);

很容易得到：

b+r=\begin{bmatrix} f(x_0)h_1/2\\ f(x_1)(h_1+h_2)/2\\ f(x_2)(h_2+h_3)/2\\ \vdots\\ f(x_{n-1}(h_{n-1}+h_n))/2\\ f(x_n)h_n/2 \end{bmatrix}+ \begin{bmatrix} \kappa_0g_0\\ \\ \\ \vdots\\ \\ \kappa_1g_1 \end{bmatrix}

MATLAB代码如下：

function b = LoadVec1D(x,f,kappa,g)
n = length(x)-1;
b = zeros(n+1,1);
for i = 1:n
    h = x(i+1) - x(i);
    b(i) = b(i) + f(x(i))*h/2;
    b(i+1) = b(i+1) + f(x(i+1))*h/2;
b(1) = b(1) + kappa(1)*g(1);
b(n+1) = b(n+1) + kappa(2)*g(2);

举个例子：假设我们现在有一维导热方程如下：

\frac{\partial ^2T}{\partial x^2}=0.03(x-6)^4, \quad x\in I=(2,8)

边界条件为：

T\big|_{x=2}=-1

\frac{\partial T}{\partial x}\big|_{x=8}=0

也就是左端是恒定温度边界条件，右端是绝热边界条件。对照之前我们理论推导，设定的边界条件应该具有下列形式：

\frac{\partial T}{\partial x}\big|_{x=2}=\kappa_0(T|_{x=0}-g_0)

\frac{\partial T}{\partial x}\big|_{x=8}=\kappa_0(T|_{x=8}-g_1)

两者貌似不一致，怎么回事？——我们之所以用诺依曼边界条件进行推导，是因为微分方程弱形式进行分部积分的时候会用到诺依曼条件，那假如是别的边界条件怎么办呢？比如本例中的就有狄利克雷条件： T\big|_{x=2}=-1 ，这就需要一点数学技巧了！

我们知道，物理模型的量一般都是有界的，不会漫天飞，因此，我们可以认为 \frac{\partial T}{\partial x}\big|_{x=2} 也是一个有限值，现在假如 \kappa_0 非常大，比如说达到 10^6 ，那 T|_{x=0}-g_0 要非常接近0才行，于是 T|_{x=0}\simeq g_0 ，所以只要我们设置 g_0=-1 即可；

第二个边界条件本身就是诺依曼条件，就简单了，直接设置 \kappa_1=0 即可；

下面我们通过MATLAB来看一下，热源的代码为：

function y = Source(x)
y = 0.03*(x-6)ˆ4; % heat source

求解结果如下：

%%solve Poission equations 
h = 0.1; % mesh size
x = 2:h:8; % mesh