20220709修改：后续的悬臂梁以及均布载荷的求算已经加上啦！同时我将它做成了matlab的APP，详情转到我的新博客：

【材料力学】基于Matlab APP Designer 开发的绘制静定梁剪力、弯矩与挠曲线的软件

梁的弯曲变形时材料力学中十分重要的研究内容，我自行设计了matlab代码以求算通过铰支座固定的梁，在多个外加集中力或是集中力偶矩的作用下，绘制梁统一的挠曲线（悬臂梁与集中载荷功能后续开发）。如有错误，欢迎交流指正。

如下示意图所示，为本程序输入参数的展示图：

相应的，程序中的输入变量如下：

程序可设定梁的弯曲刚度EI ，单位为 N/m，梁的长度单位为m，Position_Zero 为统一以梁的左端为坐标原点时两个铰支座的坐标（如上图中标出）。 Angle_Zero 为转角为零的点，如悬臂梁的末端等，因程序中暂未编写悬臂梁的求解，因此此项暂设为0。Single_Force为集中力的输入信息，矩阵第一列为集中力的施加点坐标，第二列为集中力的大小，集中力规定为向下为正。数组允许输入多行，每行相当于一个集中力。Single_Torque矩阵的输入方式类似，将集中力偶矩的正方向规定为顺时针。完成输入后，运行程序，会解算出当前载荷情况下梁的弯矩图与相应的挠曲线。

程序采用了微元化的思想，将梁的总长度划分为50000份微元（调节程序中的fineness参数来变换微元数量，当梁较长时，可以根据需要提高微元数目）。而使用者输入的外载荷的位置可能不在微元划分后的“标准坐标点”上，因此，有必要使程序有能力自主修复使用者的输出参数。为了达到这一目的，同时为了便于后续计算过程的展开，在完成输入变量后，程序将得到以下预处理信息：

其中division_Value即为微元的长度，基于这些信息，设计的坐标点“标准化”函数如下：

该函数会将input矩阵的第一列标准化至标准的微元坐标上。即实现作用点的“微移动”，这样的微移动对实际的计算结果实际上没有影响。至此，进一步调用axis_regularize 函数即可完成对输入载荷的规范化：

在完成了上述的标准化之后，便可以开始求解铰支座的反力，因为铰支座的反力都是沿铅直方向，因此统一通过“对A支座取矩求单个外载荷作用下的B支座力 -> 所求的所有B支座力叠加 -> 利用 求解A支座力”的流程来实现，具体代码如下：

至此，可以将所有支座撤去，得到纯受力图。

梁的弯矩求算则采取叠加原理的思想对包括铰支座支座力在内的所有载荷逐一求算其对梁的弯矩的影响，使用统一行向量叠加计算的方式直接完成叠加计算：

上述代码执行完毕后，即可得到长度为fineness的行向量，其中的弯矩值与x行向量中的坐标值一一对应。

因为： 则对弯矩值做两次积分，除掉梁的弯曲刚度EI ，即可得到梁的挠度值，为了提高微元算法的准确度，我使用了当前微元面积与前后共三个微元面积和的平均值作为当前位置的积分微元值，两次积分操作如下：

但仅仅进行两次积分是不够的，这相当于求解了两次数值不定积分，忽视了积分过程中产生的常数项，即真正的挠度应表示为： 而两个待定系数a、b则通过带入铰支座处的坐标求得（铰支座的挠度为0）：

至此，便完成了整个挠度的计算过程。

这部分就是较为简单的绘图功能了，不做详细介绍，具体代码如下：