声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4282|回复: 24

[编程技巧] 求助---高手帮忙调一下程序,主要求解特征值

[复制链接]
发表于 2014-7-4 17:28 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
    编了个调用程序,当I=3,J=1,3,5时,所得结果前几个数据与已知结果完全相同,当I=3,J=7,9,11时,前几个数据都有误差,请高手帮忙看看,是不是程序某个地方没有编好,可以适当调整一下吗?谢谢,文件已上传

主程序.rar

1.18 KB, 下载次数: 5

主程序

子程序.rar

2.17 KB, 下载次数: 4

调用程序

回复
分享到:

使用道具 举报

 楼主| 发表于 2014-7-4 17:57 | 显示全部楼层
下面结果都是已知结果:
I=3,J=1,得到 0.3897、1.687、2.240、3.506、4.371;
I=3,J=3,得到 0.3142、1.393、1.825、2.449、3.144;
I=3,J=5,得到 0.3100、1.374、1.810、2.197、3.124;
程序计算结果与上面的结果一一对应,而且数值一样。

程序计算结果与下面的结果不再对应,但是只要可以得到下面的数据就为正确:
根据程序所得结果,只能与下面结果接近,不能完全一致,请帮忙调一下是不是有问题:

I=3,J=7,得到  0.3092、1.370、1.806、2.167、3.113;
I=3,J=9,得到  0.3090、1.369、1.806、2.158、3.106;
I=3,J=11,得到 0.3089、1.369、1.805、2.155、3.103.

非常感谢
发表于 2014-7-11 00:06 | 显示全部楼层
看LZ连发了几个帖子, 看来好像应该都是这个问题吧!?
我是下载过了,个人水平有限但还是说说个人看法供参考
1.附件中程序的註解,我这边打开都成乱码了,无法理解学习。不长,建议考量直接贴即可, 别人也不必花体能下载
2.建议考量附上简单推导过程,方便别人理解LZ的程序
3.儘可能简化说明下整体流程

我想LZ在整理过程或许也能自我发现问题, 不然我想可能比较没人会花太多时间去关注
总归整体只有LZ最熟, 想办法让愿意帮忙的人花最少时间吧
 楼主| 发表于 2014-7-11 11:37 | 显示全部楼层
主程序如下:

clear all
format short g
syms x y
I=3;
J=7;
n=2;
L=0.5;     
R=1;      
Ri=0.2;
H=0.2;     
mu=0.3;  
t=2*mu/(1-2*mu);  

T1=pi;
T2=pi;

%%%%无量纲值计算公式
r=Ri/R;
h=H/R;
l=R/L;

yi=@(x) l*sqrt((1-h/2)^2-x.^2);
yo=@(x) l*sqrt((1+h/2)^2-x.^2);

%%%%质量、刚度子矩阵
M=zeros(3*(I+1)*(J+1));
m11=zeros((I+1)*(J+1));
m22=m11;
m33=m11;
m12=m11;
m13=m11;
m21=m11;
m23=m11;
m31=m11;
m32=m11;

K=zeros(3*(I+1)*(J+1));
k11=zeros((I+1)*(J+1));
k12=k11;
k13=k11;
k21=k11;
k22=k11;
k23=k11;
k31=k11;
k32=k11;
k33=k11;

%%%%求解质量子矩阵、刚度子矩阵内的各个元素
for  iii=1:(I+1)*(J+1)
    for  jjj=1:(I+1)*(J+1)
   
          jj=mod(jjj,J+1)-1;     %取模(余)
       if  jj==-1
           jj=J;
       end

          j=mod(iii,J+1)-1;
       if  j==-1
           j=J;
       end

          ii=floor((jjj-jj)/(J+1));   %取整
         
           i=floor((iii-j)/(J+1));
        
%%%%调用各个函数
f1=f_1(i,j,ii,jj,yi,yo,r,h);
f2=f_2(i,j,ii,jj,yi,yo,r,h);
f3=f_3(i,j,ii,jj,yi,yo,r,h);
f4=f_4(i,j,ii,jj,yi,yo,r,h);
f5=f_5(i,j,ii,jj,yi,yo,r,h);
f6=f_6(i,j,ii,jj,yi,yo,r,h);
f7=f_7(i,j,ii,jj,yi,yo,r,h);
f8=f_8(i,j,ii,jj,yi,yo,r,h);
f9=f_9(i,j,ii,jj,yi,yo,r,h);
f10=f_10(i,j,ii,jj,yi,yo,r,h);

%%%%%计算质量子矩阵、刚度子矩阵内的各个元素
m11(iii,jjj)=T1*f10;
m22(iii,jjj)=T1*f10;
m33(iii,jjj)=T2*f10;

k11(iii,jjj)=T1*((t+2)*(f1+f2)+t*(f4+f5)+l^2*f3)+n^2*T2*f1;
k12(iii,jjj)=l*T1*(t*(f7+f9)+f8);
k13(iii,jjj)=n*T1*((t+2)*f1+t*f4)+n*T2*(f1-f5);

%k21(iii,jjj)=l*T1*(t*(f6+f8)+f9);
k22(iii,jjj)=T1*((t+2)*l^2*f3+f2)+n^2*T2*f1;
k23(iii,jjj)=n*l*(t*T1*f6-T2*f7);

%k31(iii,jjj)=n*T1*((t+2)*f1+t*f5)+n*T2*(f1-f4);
%k32(iii,jjj)=n*l*(t*T1*f7-T2*f6);
k33(iii,jjj)=n^2*(t+2)*T1*f1+T2*(l^2*f3+f1+f2-f4-f5);

    end
end
k21=k12';
k31=k13';
k32=k23';

%%形成总质量、总刚度矩阵
M=[m11,m12,m13;m21,m22,m23;m31,m32,m33];
K=[k11,k12,k13;k21,k22,k23;k31,k32,k33];
MM=balance(M);  %相似变换,特征值相同,精度更高
KK=balance(K);
[V1,Ommiga1]=eig(KK,MM,'qz');  
%[V1,Ommiga1]=eig(K,M,'qz');  
Ommiga2=diag(Ommiga1);
Ommiga3=sqrt(Ommiga2);         
for m=1:length(Ommiga3)
    w(m)=isreal(Ommiga3(m));  %实数特征值所在的位置
end
V2=V1(:,w);       %过滤虚数特征值后所得到的特征向量
[Ommiga3,n]=sort(Ommiga3(w));
V3=V2(:,n);       %实数特征值排序后对应的特征向量
ommiga=Ommiga3
 楼主| 发表于 2014-7-11 11:39 | 显示全部楼层
各个子程序如下:
function z=f_1(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) x.^(a+c-1).*y.^(b+d),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) x.^(a+c-1).*y.^(b+d),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_2(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) a*c*x.^(a+c-1).*y.^(b+d),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) a*c*x.^(a+c-1).*y.^(b+d),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_3(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) b*d*x.^(a+c+1).*y.^(b+d-2),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) b*d*x.^(a+c+1).*y.^(b+d-2),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_4(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) a*x.^(a+c-1).*y.^(b+d),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) a*x.^(a+c-1).*y.^(b+d),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_5(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) c*x.^(a+c-1).*y.^(b+d),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) c*x.^(a+c-1).*y.^(b+d),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_6(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) b*x.^(a+c).*y.^(b+d-1),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) b*x.^(a+c).*y.^(b+d-1),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_7(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) d*x.^(a+c).*y.^(b+d-1),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) d*x.^(a+c).*y.^(b+d-1),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_8(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) c*b*x.^(a+c).*y.^(b+d-1),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) c*b*x.^(a+c).*y.^(b+d-1),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_9(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) a*d*x.^(a+c).*y.^(b+d-1),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) a*d*x.^(a+c).*y.^(b+d-1),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

function z=f_10(a,b,c,d,yi,yo,r,h)
z1=quad2d(@(x,y) x.^(a+c+1).*y.^(b+d),r,1+h/2,-1,yo);
z2=quad2d(@(x,y) x.^(a+c+1).*y.^(b+d),r,1-h/2,-1,yi);
z=z1-z2;
z=vpa(z);

发表于 2014-7-11 12:26 | 显示全部楼层
代码看起来有点费劲。

对于原理和方法,可以先整理一下,放上来看看。

点评

赞成: 5.0
赞成: 5
  发表于 2014-7-11 13:57
发表于 2014-7-11 12:52 | 显示全部楼层
其实你算出来的结果中,还不止你列出的这几个,中间还有一些结果,可能是约束问题。
 楼主| 发表于 2014-7-11 15:53 | 显示全部楼层
这个是自由约束的结果
 楼主| 发表于 2014-7-11 16:03 | 显示全部楼层
这两个图是结构的三个方向位移假定,r、z和θ是径向、轴向和环向三个方向坐标(圆柱坐标系),I=K=M,J=L=N。
~`8N@L{ZLWZH%6RY2B)HQ1P.jpg
4QINW$)S~7%V[UBJSDRI2@V.jpg
 楼主| 发表于 2014-7-11 16:12 | 显示全部楼层
原理就是通过瑞雷-里兹法 假定结构的位移函数(如上图),然后求出结构的最大动能和势能,利用能量变分得到刚度矩阵、质量矩阵的各个元素,只是此时每个元素都是个子矩阵,多项式项数越多,次数也越多,得到的矩阵阶数都在增加,当I=3、J=7的时候,得到的总刚度、质量矩阵行数 列数都是3*(3+1)*(7+1)=96的矩阵,然后求解频率特征值
 楼主| 发表于 2014-7-11 16:23 | 显示全部楼层
dollfish000 发表于 2014-7-11 12:52
其实你算出来的结果中,还不止你列出的这几个,中间还有一些结果,可能是约束问题。

因为里兹法求解出来的一系列频率,不一定跟结构的真实频率的每一阶一一对应的,只有位移项数取得足够多的时候,也许可以包括结构的所有阶频率值
发表于 2014-7-11 16:36 | 显示全部楼层
从开始假设这个模态,到最后计算并不能包含左右的频率成份,连最低阶的几个都不能包含,我不觉得有什么意义;另外像这个模型,做结构工程的人一看到这个就觉得非常简单,而且你说是算自由振动没有任何约束,那更加好求了,就是一根梁。程序上我也看不出有什么问题,毕竟编程的风格不太一样。
 楼主| 发表于 2014-7-11 16:49 | 显示全部楼层
dollfish000 发表于 2014-7-11 16:36
从开始假设这个模态,到最后计算并不能包含左右的频率成份,连最低阶的几个都不能包含,我不觉得有什么意义 ...

这是个圆柱壳结构,频率跟环向半波数n有关系,不同的n,都可以得到一系列频率,从低到高,位移函数不是很先进,看似简单,主要就是推导过程比较新(得到动能和势能的方法很好),从而得到的矩阵各个元素也不是很难
 楼主| 发表于 2014-7-11 16:50 | 显示全部楼层
由于得到的刚度和质量矩阵的元素下标都是4个字母的,(MATLAB中矩阵元素只有两个),所以先找到四个下标i、j、ii及jj同两个下标iii和jjj的关系,如图,M和K矩阵的每个元素都是下标i、j、ii及jj的函数,即f1-f10十个函数。所以每次得到一个元素都要调用一次函数(10个),所以每个元素就是一个子矩阵
0PEKC8CQ4@60SH@`T%7`YNL.jpg
发表于 2014-7-11 22:18 | 显示全部楼层
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-25 01:06 , Processed in 0.062661 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表