本文介绍分别用Matlab，COMSOL， 以及自己手工 求解拉普拉斯方程。
1.Matlab
命令行输入
1. pdetool
2. 启动界面，工具栏选择画长方形工具，画一个正方形，如图1-1
3. 点击菜单 Options--》Application--》Heat transfer
4. 点击菜单 Boundary--》 Boundary Mode， 选中图中左边一条边，温度设置为100， 选中右边一条边温度设为0.
5. 点击菜单 Solve--Solve PDE 如图1-2
图1-1：
图1-2：
Matlab 解偏微分方程最大的局限性是 只能求解平面不能解三维。
2. COMSOL
1. Space 选择2D
2. 求解类型选择Heat transfer-》Heat transfer in solid，求解类型选Stationary
3. 选中树形菜单Geometry，右键，选择rectangle，然后画出一个正方形
4. 选中树形菜单Material，选中一个Build-in的材料
5. 选中树形菜单Heat Transfer in Solids,右键选中 Temperature，然后选中长方形左边设置100,
6. 重复5，选中正方形右边设置0
7. 点击树形菜单Study1，右键计算 Compute selected step
结果看起来和Matlab 不太一样，是因为材料属性热导率不同，Matlab 默认边界温度为0。
3. 手工求解
先将正方形离散为三角单元，为三角单元建立形函数，然后利用伽辽金方法推导出积分公式，加入边界，建立可编程实现的有限元方程，求解线性方程，计算出结果。
三角单元离散化参考：FEM之单元(1)---三角单元介绍 。
假设三角单元温度场函数为：
T = [S1 S2 S3] [Ti Tj Tk]；
S1 S2 S3为 线性形函数；
S1=(a1+b1x+c1y)/2A；
S2=(a2+b2x+c2y)/2A；
S3=(a3+b3x+c3y)/2A；
Ti Tj Tk 为 三角形三个顶点上的温度；
A 为三角形面积；
T 表示了三角单元内部任意一点的温度。
下面最主要的一步是利用伽辽金方法为三角单元建立方程：
将形函数与偏微分方程乘积的积分余差为0.
其中
即为[S1 S2 S3]形函数，其中包含了变量x，y，上式可以拆成两个
对第一个式子进行解析，按照分步积分和链式求导法则，有如下表达式：
将该式子右边第二项移动到左边，就得到了要计算的表达式其中
求导无任何问题，简单的高等数学知识；表达式
计算要利用一下 格林函数，将面积分转化成曲线积分，其实也是大学高等数学知识。最终的计算结果可以将积分表达式转化成用a，b，c表现的矩阵形式。
上面一步是很多教科书都没有写明的，一带而过，但这才是用伽辽金方法生成有限元方程的关键所在。
形成矩阵形式以后，就可以很方便的用编程实现了。详细推导过程可以参考《有限元分析--ANSYS理论与应用》P310--P312， 作者 Saeed Moavenl，这个系列的书都不错，深入浅出，通俗易懂。
下一步就是要将单个矩阵组装成整体矩阵表达式。知道了节点编号，就是对号入座就行了。假设节点数目为N，生成一个N*N的总体矩阵K，按节点编号，将单个单元标号的矩阵放入整体矩阵中，不同单元相同节点号直接相加即可，这表明了单元的连续性。之后就是边界条件，将已知边界的温度生成一个N的向量B，解方程 Kx=B，即可得到温度场分布。在解决实际问题时，由于矩阵特别大，在总刚组装以及求解线性方程组时会非常讲究，以后有时间讨论。
注：文中所使用的软件仅供学习使用。
