24小时热门版块排行榜

返回列表

spvfly

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1131
帖子: 148
在线: 54.3小时
虫号: 812244
注册: 2009-07-20
性别: GG
专业: 应用地球物理学

[求助] 一段用MATLAB编写的代码，结果错误，请各位帮忙找找错误出在哪了！！！已有1人参与

%%%主程序
XY=xlsread('data.xlsx','ND');
I3=xlsread('data.xlsx','NE');
UB=xlsread('data.xlsx','UB');
ND=33;%点数
NE=40;%单元数
ND1=24;%边值点数
UK=UK1(ND,NE,I3,XY);%总系数矩阵
[UK,UB1]=UBK(UK,UB,ND1);%将KU=0方程组变为KU=B形式的方程组
U=UK\UB1'%计算得到结果
%%%%%%%%
总系数矩阵程序
%%%%%%%%
function y=UK1(ND,NE,I3,XY)
X=[0,0,0];
Y=[0,0,0];
SK=zeros(ND,ND);
KE=zeros(3,3);
for P=1:NE %遍历各个单元
      for I=1:3 %遍历一个单元的三个节点
         L=I3(I,P); %确定一个单元的节点
         X(I)=XY(1,L); %该节点的x坐标
         Y(I)=XY(2,L); %该节点的y左边
      end
      KE=UKE1(X,Y);%一个单元的单元系数矩阵
      %下面这段程序是为了把单元系数矩阵放入定带宽矩阵中去
      for J=1:3
         for K=1:3
            NJ=I3(J,P);
            NK=I3(K,P);
            SK(NJ,NK)=SK(NJ,NK)+KE(J,K);
         end
      end
end
y=SK;
%%%%%%%
单元系数矩阵
%%%%%%%
function y=UKE1(X,Y)
K1=zeros(3,3);%定义一个3*3数组用于存放单元系数矩阵的值
A=[0,0,0];
B=[0,0,0];
C=[0,0,0];
A(1)=Y(2)-Y(3);
A(2)=Y(3)-Y(1);
A(3)=Y(1)-Y(2);
B(1)=X(3)-X(2);
B(2)=X(1)-X(3);
B(3)=X(2)-X(1);
S=2*(A(1)*B(2)-A(2)*B(1));
for I=1:3
for J=1:3
      K1(I,J)=(A(I)*A(J)+B(I)*B(J))/S;
end
end
y=K1;
%%%%%%%%%
将KU=0形式转换为KU=B形式
%%%%%%%%%
function [X,Y]=UBK(UK,UB,ND1)
UB1=zeros(1,33);
for I=1:ND1
D=UB(1,I);
V=UB(2,I);
UK(D,D)=UK(D,D)*1e10;
UB1(D)=V*UK(D,D);
end
X=UK;
Y=UB1;

%%%%%%%
所得到的错误结果
%%%%%%
程序结果解析解
0.354616826 0.302042805
0.455989065 0.355568018
0.562915374 0.41747372
0.675049586 0.482845483
0.767370579 0.526986487
0.761868629 0.517705616
0.621009799 0.455065182
0.46388813 0.373490832
0.342221541 0.301003461

回复此楼