| 查看: 750 | 回复: 3 | |||
[求助]
一段用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 |
回复此楼» 猜你喜欢
散金币祈福
已经有186人回复
-
求助，如何提取ELK的rt-TDDFT在某一时刻的自旋密度分布和电子密度分布
已经有1人回复
-
物理学I论文润色/翻译怎么收费?
已经有240人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
澳门大学IAPME李海峰教授常年招收博士研究生
已经有0人回复
-
深圳大学应用技术学院招聘凝聚态物理博士后
已经有112人回复
-
深圳大学2027级光电信息工程/物理学研究生推免
已经有69人回复
-
求助中医院临床中医中标的青年标书学习，谢谢
已经有4人回复
-
浙大物院表面物理课题组长期招收普博及直博学生
已经有257人回复
-
【征稿】Springer Nature热电专题 Discover Materials / Discover Chemistry
已经有39人回复
2楼2016-12-25 22:10:38
3楼2016-12-26 20:48:25
FMStation
至尊木虫 (知名作家)
- 应助: 591 (博士)
- 贵宾: 0.03
- 金币: 18761.2
- 红花: 97
- 帖子: 8891
- 在线: 1497.8小时
- 虫号: 2400059
- 注册: 2013-04-04
- 专业: 计算机应用技术
4楼2016-12-28 20:38:37
20