求助---高手帮忙调一下程序，主要求解特征值

江山美人

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

江山美人

下面结果都是已知结果：
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.

非常感谢

ChaChing

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

我想LZ在整理过程或许也能自我发现问题, 不然我想可能比较没人会花太多时间去关注
总归整体只有LZ最熟, 想办法让愿意帮忙的人花最少时间吧

江山美人

主程序如下：

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

江山美人

各个子程序如下：
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);

wanyeqing2003

代码看起来有点费劲。

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

dollfish000

其实你算出来的结果中，还不止你列出的这几个，中间还有一些结果，可能是约束问题。

江山美人

这个是自由约束的结果

江山美人

这两个图是结构的三个方向位移假定，r、z和θ是径向、轴向和环向三个方向坐标（圆柱坐标系），I=K=M，J=L=N。

江山美人

原理就是通过瑞雷-里兹法 假定结构的位移函数（如上图），然后求出结构的最大动能和势能，利用能量变分得到刚度矩阵、质量矩阵的各个元素，只是此时每个元素都是个子矩阵，多项式项数越多，次数也越多，得到的矩阵阶数都在增加，当I=3、J=7的时候，得到的总刚度、质量矩阵行数 列数都是3*（3+1）*（7+1）=96的矩阵，然后求解频率特征值

江山美人

dollfish000 发表于 2014-7-11 12:52
其实你算出来的结果中，还不止你列出的这几个，中间还有一些结果，可能是约束问题。

因为里兹法求解出来的一系列频率，不一定跟结构的真实频率的每一阶一一对应的，只有位移项数取得足够多的时候，也许可以包括结构的所有阶频率值

dollfish000

从开始假设这个模态，到最后计算并不能包含左右的频率成份，连最低阶的几个都不能包含，我不觉得有什么意义；另外像这个模型，做结构工程的人一看到这个就觉得非常简单，而且你说是算自由振动没有任何约束，那更加好求了，就是一根梁。程序上我也看不出有什么问题，毕竟编程的风格不太一样。

江山美人

dollfish000 发表于 2014-7-11 16:36
从开始假设这个模态，到最后计算并不能包含左右的频率成份，连最低阶的几个都不能包含，我不觉得有什么意义 ...

这是个圆柱壳结构，频率跟环向半波数n有关系，不同的n，都可以得到一系列频率，从低到高，位移函数不是很先进，看似简单，主要就是推导过程比较新（得到动能和势能的方法很好），从而得到的矩阵各个元素也不是很难

江山美人

由于得到的刚度和质量矩阵的元素下标都是4个字母的，（MATLAB中矩阵元素只有两个），所以先找到四个下标i、j、ii及jj同两个下标iii和jjj的关系，如图，M和K矩阵的每个元素都是下标i、j、ii及jj的函数，即f1-f10十个函数。所以每次得到一个元素都要调用一次函数（10个），所以每个元素就是一个子矩阵

qjf880726