这是第一个运算过程：
clc,clear;
E=2e11;A=6.65e-2;I=118.6e-6;q=25000;
%输入弹性模量，单位为（pa）；横截面积，单位为（m^2）；极惯性矩，单位为(m^4)；等效载荷，单位为（N/m）
O=[1 2;2 3;3 4;4 5;5 6;6 7];
%输入单元编号与单元节点编号矩阵
V=[1 2;3 4;5 6; 7 8;9 10;11 12;13 14];
%输入节点自由度编号矩阵
B=[0 0;1250 0;2500 0;3750 0;5000 0;6250 0;7500 0];
%输入节点坐标矩阵，其单位为（mm）
K=zeros(4,4);ZK=zeros(14,14);T=zeros(14,14);
%定义三个要求或之后要应用的矩阵
u=[0 0 1 1 1 1 1 1 1 0 1 1 1 1]';
%定义位移矩阵，未知项令其为1，方便之后计算
F1=[-0.5q1.25 -q*(1.25^2)/12 -q1.25 0 -q1.25 0 -q1.25 0 -q1.25 0 -q1.25 0 -0.5q1.25 q(1.25^2)/12]';
%定义外力矩阵，其单位为（N）
for n1=1:6
G(1:1,:)=B(O(n1,1):O(n1,1),:);
G(2:2,:)=B(O(n1,2):O(n1,2),:);
%定义一个过渡矩阵用于传入单元刚度矩阵子函数，便于计算与传入
K=Untitled12(G,E,I);
%调用单元刚度矩阵计算函数
DOF(1)=V(O(n1,1),1);
DOF(2)=V(O(n1,1),2);
DOF(3)=V(O(n1,2),1);
DOF(4)=V(O(n1,2),2);
%定义单元坐标与总刚度矩阵坐标转化的DOF矩阵，便于之后总体刚度矩阵组集
ZK=Untitled13(K,DOF,ZK);
%调用总体刚度矩阵组集函数
end
Q=Untitled14(ZK,u,F1);
%调用位移计算函数计算位移，其单位为（m）
P=ZK*Q;
%通过总体刚度矩阵于位移矩阵计算出等效节点载荷矩阵，其单位为（N）
F=P-F1;
%通过等效节点载荷与外力矩阵计算支反力，其单位为（N）
for n3=1:6
G1(1:1,:)=B(O(n3,1):O(n3,1),:);
G1(2:2,:)=B(O(n3,2):O(n3,2),:);
%过渡矩阵，与之前作用相同
Q1=[Q(V(O(n3,1),1));Q(V(O(n3,1),2));Q(V(O(n3,2),1));Q(V(O(n3,2),2))];
SS=Untitled15(E,Q1,G1);
%通过应力应变关系求出应力，其单位为（Pa）
SK(n3)=SS;
%SK矩阵为储存应力的矩阵
end
function [K1]=Untitled12(G,E,I)
x1=G(1,1);y1=G(1,2);x2=G(2,1);y2=G(2,2);
%设置通用节点坐标
L=sqrt((x2-x1)^2+(y1-y2)^2);
%通过坐标计算单元杆长度
K1=EI/(L^3)[12 6L -12 6L;6L 4L^2 -6L 2L^2;-12 -6L 12 -6L;6L 2L^2 -6L 4L^2];
%单元刚度矩阵计算
end
function [T]=Untitled13(K,DOF,ZK)
for n1=1:4
for n2=1:4
ZK(DOF(n1),DOF(n2))=ZK(DOF(n1),DOF(n2))+K(n1,n2);
end
end
T=ZK;
end
function [U]=Untitled14(ZK,u,F1)
%使用置一法计算位移边界条件
for n1=1:14
if(u(n1,1)==0)%判断位移是否为0
for n2=1:14
ZK(n1,n2)=0;
end
for n3=1:14
ZK(n3,n1)=0;
end
ZK(n1,n1)=1;
F1(n1,1)=0;
%运用置一法编辑整体刚度矩阵
end
end
U=inv(ZK)*F1;
%通过总体刚度矩阵逆矩阵与外力矩阵计算位移矩阵
end
function [stress]=Untitled15(E,Q1,G1)
x1=G1(1,1);x2=G1(2,1);
%设置通用节点坐标
x=(x1+x2)/2;y=0.1585;
%设所求点为单元上表面中点
L=1.25;
%输入单元长度
e=x/L;
%输入所求点与单元左节点的距离与单元长度的比例
B1=(12e-6)/(LL);
B2=(6e-4)/L;
B3=-(12e-6)/(LL);
B4=(6e-2)/L;
B=-y*[B1,B2,B3,B4];
%输入几何矩阵
stress=EBQ1;
%计算该单元在所求点的应力
end
这是第二个运算过程：
clc;
clear;
E=2e11;%输入弹性模量
A=6.65e-2;%横截面积
I=118.6e-6;%极惯性矩
q=25000;%荷载
bh=[1 2;2 3;3 4;4 5;5 6;6 7];%单元编号
zb=[0 0;1250 0;2500 0;3750 0;5000 0;6250 0;7500 0];%节点坐标矩阵
zk=zeros(14,14);
u=[0 0 1 1 1 1 1 1 1 0 1 1 1 1]';
f1=[-0.5q1.25 -q*(1.25^2)/12 -q1.25 0 -q1.25 0 -q1.25 0 -q1.25 0 -q1.25 0 -0.5q1.25 q(1.25^2)/12]';
for n1=1:6
gd(1,:)=zb(bh(n1,1),:);
gd(2,:)=zb(bh(n1,2),:);
Ke=danyuan(gd,E,I);
zk=zuji(Ke,zk,n1,bh);
end
r=weiyi(zk,u,f1);
%调用位移计算函数计算位移，其单位为（m）
P=zkr;
%通过总体刚度矩阵于位移矩阵计算出等效节点载荷矩阵，其单位为（N）
F=P-f1;
%通过等效节点载荷与外力矩阵计算支反力，其单位为（N）
for n2=1:6
gd1(1,:)=zb(bh(n2,1),:);
gd1(2,:)=zb(bh(n2,2),:);
Q=[r(2bh(n2,1)-1);r(2bh(n2,1));r(2bh(n2,2)-1);r(2*bh(n2,2))];
stress=yingli(E,Q,gd1);%一次循环得到一个
sk(n2)=stress;%总循环完毕得到6个组成1x6的矩阵
end
function [ Ke ] = danyuan( gd,E,I )
%输入坐标过渡矩阵gd，弹性模量，极惯性矩
%输出单元坐标刚度矩阵
x1=gd(1,1);y1=gd(1,2);
x2=gd(2,1);y2=gd(2,2);
L=sqrt((x2-x1)^2+(y2-y1)^2);
Ke=EI/(L^3)[12 6L -12 6L;6L 4L^2 -6L 2L^2;-12 6L 12 -6L;6L 2L^2 -6L 4L^2];
end
function[zk]=zuji(Ke,zk,n1,bh)
%该函数进行单元刚度矩阵的组装
%输入清零的结构总的刚度矩阵zk,自由度编号DOF,单元刚度矩阵Ke
%输出组集好的总刚度矩阵zk
L=length(Ke);%L为zk1矩阵行数或者列数中的较大值
DOF(1)=2bh(n1,1)-1;
DOF(2)=2bh(n1,1);
DOF(3)=2bh(n1,2)-1;
DOF(4)=2bh(n1,2);
for a=1:L
for b=1:L
zk(DOF(a),DOF(b))=zk(DOF(a),DOF(b))+Ke(a,b);%将单元刚度矩阵Ke的元素依次带入总体刚度矩阵zk中
end
end
function[q1]=weiyi(zk,u,f1)
%该函数计算结构的节点位移和节点总载荷矩阵
%输入整体的刚度矩阵zk，节点的外力列阵f1，结构的节点外载列阵p
%输出结构的节点位移q1
for n1=1:14
if(u(n1,1)==0)%判断位移是否为0
for n2=1:14
zk(n1,n2)=0;
end
for n3=1:14
zk(n3,n1)=0;
end
zk(n1,n1)=1;
f1(n1,1)=0;
%运用置一法编辑整体刚度矩阵
end
end
q1=inv(zk)*f1;
end
function [stress]=yingli(E,Q,gd1)
%该函数计算单元的应力
%输入弹性模量E，自由度编号Q,单个单元坐标gd1
%输出单元应力
x1=gd1(1,1);x2=gd1(2,1);
x=(x1+x2)/2;y=0.1585;
L=1.25;
e=x/L;
B1=(12e-6)/(LL);
B2=(6e-4)/L;
B3=-(12e-6)/(LL);
B4=(6e-2)/L;
B=-y*[B1 B2 B3 B4];
stress=EBQ;
end
两个计算获得的zk以及f1完全相同，但在inv(zk)*f1得出的结果却不一致，求大佬解答
请把代码按照代码格式粘贴（参见编辑框右上角工具栏</>）