MATLAB报告
Matlab程序求解简要过程如下：
（1）求取单元节点位移提取矩阵T
单元节点位移提取矩阵T本质上是置换矩阵群中的一个，结果可将任意杂乱的节点顺序置换成统一的顺序。

另一方面其作用是对单元刚度矩阵进行“升维操作”，将单元刚度矩阵统筹到整体刚度矩阵上来，便于对总体节点位移矩阵和支座反力进行求取。

本程序分析过程中对单元1的节点提取是按顺序编号1-2-3，对单元2的节点提取是按顺序编号2-3-4。

单元1的节点位移提取矩阵如下：
单元2的节点位移提取矩阵如下：
（2）求取单元几何矩阵B
单元1的节点按编号顺序1-2-3分别进行对几何函数矩阵或算子矩阵的bi逆时针操作，对ci顺时针操作；单元2的节点按编号顺序2-3-4分别进行对几何函数矩阵的bi顺时针操作，对ci逆时针操作.在MATLAB程序中通过mod()取模函数来达到对节点的顺时针或逆时针循环操作。

单元1的几何矩阵如下：
单元2的几何矩阵如下：
（3）求取应力矩阵S
单元应力矩阵满足S=D*B,其中D为弹性矩阵，B为单元几何矩阵
各单元的弹性矩阵如下：
单元1的应力矩阵如下：
单元2的应力矩阵如下：
（4）求取单元刚度矩阵K
单元刚度矩阵K满足公式K=B’*D*B*t*A,其中t为平面板的厚度，A为单元面积，且单元刚度矩阵为对称矩阵。

单元1的刚度矩阵如下：
单元2的刚度矩阵如下：
（5）求取总体刚度矩阵sumKK
由上述步骤求得的单元刚度矩阵K利用单元虚功原理和刚度方程可导出K’*δ=f，其中δ为单元节点位移列阵，f为单元等效节点载荷列阵，为了能将各个单元刚度方程统一到一个整体，便需要步骤（1）的单元节点提取矩阵对单元刚度方程进行变换，将两个变换结果联立便得到总体刚度方程，其中也可得到总体刚度矩阵sumKK，且总体刚度矩阵可由sumKK=ΣT’*K*T求得。

总体刚度矩阵如下：
（6）求取总体节点位移矩阵和支座反力
利用上述步骤提到的总体刚度方程sumKK*delta=F，其中delta为总体节点位移矩阵，F为总体等效节点载荷列阵。

利用总体边界条件可将总体刚度方程分成两部分，一部分先求取总体节点位移矩阵delta，最后代入总体刚度方程求取总体等效节点载荷列阵F，由此便得到支座反力。

总体节点位移矩阵如下：
支座反力如：
（7）求取单元应力
各单元应力满足θ=S*δ，其中S为单元应力矩阵，δ为单元节点位移列阵。

单元1的应力：
单元2的应力：
下图为应变前后的各单元对比图：
（注：可编辑下载，若有不当之处，请指正，谢谢!）。