关于化简的问题，各位前辈帮帮忙！

hanxiao

程序的简要说明：
这是一个分段梁应用振型求响应的程序。t代表输入激励的频率；z是梁的长度坐标；
a就是输入的振型；s是输入的载荷谱；中间的循环是因为阶梯型梁嘛，参数是分段变化的；
yi和yr是求得的位移响应的虚部和实部，它们是频率t和坐标z的函数；语句z=38后，
是求z=38米处的位移响应的功率谱密度函数，即s38uuu和s38uu，它们是频率t的函数；
然后我想画出这个功率谱密度函数的图像，并查看它的数据结果。

跪求各位大虾帮帮忙！
不胜感激！


%-- 06-10-20下午7:21 --%
syms z t
a=[1-cos(pi*z/(2*68)),1-cos(3*pi*z/(2*68)),1-cos(5*pi*z/(2*68)),1-cos(7*pi*z/(2*68)),1-cos(9*pi*z/(2*68))];   
a1=diff(a);                                                                              
a2=diff(a1);
s=(15.7^4+4*0.63*0.63*15.7^2*t^2)*0.001574/((15.7^2-t^2)^2+4*0.63*0.63*t^2*15.7^2);              

%关于惯性矩的循环，38次
r=0;
m2=0;
k1=0;
p2=0;
for i=1:38
b0=int(a'*a,z,i-1,i);
b1=int(a1'*a1,z,i-1,i);
b2=int(a2'*a2,z,i-1,i);
bp=int(a',z,i-1,i);
gxj=55-i;
if  i>8&i<=29
midu=22994.15;
else
midu=7800;
end
r=r+gxj*midu^2*b0;
m2=m2+gxj*midu*b1;
k1=k1+gxj*b2;
p2=p2+gxj*midu^2*bp;
end
%循环结束，得到0-38米各阵
rz38=2.47*10^(-11)*r;
m2z38=6.19*m2;
k1z38=2.1*10^11*k1;
p2z38=2.47*10^(-11)*t^2*sa*p2;
%有关惯性矩的计算38-47米(9米一段)
rz9=2.55*10^(-2)*int(a'*a,z,38,47);
m2z9=8.20*10^5*int(a1'*a1,z,38,47);
k1z9=3.57*10^12*int(a2'*a2,z,38,47);
p2z9=2.55*10^(-2)*t^2*sa*int(a',z,38,47);
%有关惯性矩的计算47-54米(第一层甲板）
rz7=42.3*int(a'*a,z,47,54);
m2z7=3.34*10^7*int(a1'*a1,z,47,54);
k1z7=3.57*10^12*int(a2'*a2,z,47,54);
p2z7=42.3*t^2*sa*int(a',z,47,54);
%有关惯性矩的计算54-68米(第二层甲板）
rz14=53*int(a'*a,z,54,68);
m2z14=3.74*10^7*int(a1'*a1,z,54,68);
k1z14=3.57*10^12*int(a2'*a2,z,54,68);
p2z14=53*t^2*sa*int(a',z,54,68);
%有关惯性矩的总阵(r,m2,k1,p2)
rz=rz38+rz9+rz7+rz14;
m2z=m2z38+m2z9+m2z7+m2z14;
k1z=k1z38+k1z9+k1z7+k1z14;
p2z=p2z38+p2z9+p2z7+p2z14;

%有关截面面积的m1
m1z1=3666*int(a'*a,z,0,8);
m1z2=22994.15*int(a'*a,z,8,29);
m1z3=4446*int(a'*a,z,29,38);
m1z38=m1z1+m1z2+m1z3;                       
m1z9=4.43*10^3*int(a'*a,z,38,47);         
m1z7=1.80*10^5*int(a'*a,z,47,54);         
m1z14=2.02*10^5*int(a'*a,z,54,68);            
%有关截面面积的p1
p1z1=-3666*int(a',z,0,8);
p1z2=-22994.15*int(a',z,8,29);
p1z3=-4446*int(a',z,29,38);
p1z38=(p1z1+p1z2+p1z3)*sa;                     
p1z9=-4.43*10^3*sa*int(a',z,38,47);              
p1z7=-1.80*10^5*sa*int(a',z,47,54);                    
p1z14=-2.02*10^5*sa*int(a',z,54,68);                     
%总阵
m1z=m1z38+m1z9+m1z7+m1z14;
p1z=p1z38+p1z9+p1z7+p1z14;

%轴力的计算0-38米(k2)
k2z=-39.15*10^6*int(a1'*a1,z,0,38);

%各阵的合并
m=simplify(m1z+m2z);
k=simplify(k1z+k2z);
p=simplify(p1z+p2z);
r=simplify(rz);
%求解计算
e1=t^4*r+k-t^2*m;
d1=-0.01*t*k;
e2=simplify(e1);
d2=simplify(d1);
e=vpa(e2,8);
d=vpa(d2,8);
i=inv(e);
j=inv(d);
q1=i*d+j*e;
q2=simplify(q1);
q3=vpa(q2,8);
q=inv(q3);
yiii=q*i*p;
yrrr=q*j*p;
yii=simplify(yiii)
yrr=simplify(yrrr)
yi=vpa(yii,8);
yr=vpa(yrr,8);


z=38;
a38=[1-cos(pi*z/(2*68)),1-cos(3*pi*z/(2*68)),1-cos(5*pi*z/(2*68)),1-cos(7*pi*z/(2*68)),1-cos(9*pi*z/(2*68))];     
a=vpa(a38,8);
urrr38=a(1)*yr(1)+a(2)*yr(2)+a(3)*yr(3)+a(4)*yr(4)+a(5)*yr(5);
uiii38=a(1)*yi(1)+a(2)*yi(2)+a(3)*yi(3)+a(4)*yi(4)+a(5)*yi(5);
urr38=simplify(urrr38)
uii38=simplify(uiii38)
ur38=vpa(urr38,8);
ui38=vpa(uii38,8);
s38uuu=ur38^2+ui38^2;
s38uu=simplify(s38uuu)
s38uuu1=2*s38uu;
s38uuu2=t^2*s38uuu1;
s38uu1=simplify(s38uuu1)
s38uu2=simplify(s38uuu2)
s3801=vpa(s38uu1,2)
s3802=vpa(s38uu2,2)

以下是我的问题：
（1)程序不复杂,计算结果十分复杂.
s38uu,s3801和s3802是我想要的计算结果,长得吓人!我想怎么能化简一下?
我是临时抱佛脚,刚刚学的matlab.想问一下 :
(2)对s38uu怎么化出它 的图形,并显示图形的数据结果?
(3)怎么对s3801和s3802进行数值积分,怎么查看积分结果?

[ 本帖最后由 hanxiao 于 2006-12-19 10:13 编辑 ]

eight

原帖由 hanxiao 于 2006-12-7 10:59 发表
程序的简要说明：
这是一个分段梁应用振型求响应的程序。t代表输入激励的频率；z是梁的长度坐标；
a就是输入的振型；s是输入的载荷谱；中间的循环是因为阶梯型梁嘛，参数是分段变化的；
yi和yr是求得的位移响应 ...

建议把题目修改，不然要不被封，要不帖子被删除。当然了，你可以选择把跪的图片附上

lxq

以后最好不用这样的标题：（

hanxiao

嗬嗬嗬嗬，知道了：）
主要是比较着急么！

xjzuo

请讲清楚一下你要实现的计算.
比如下面这一段代码你想算什么(b0,b1,...是想存储那些部分等):
%%%%%%%%%%%%%%
for i=1:38
b0=int(a'*a,z,i-1,i);
b1=int(a1'*a1,z,i-1,i);
b2=int(a2'*a2,z,i-1,i);
bp=int(a',z,i-1,i);
gxj=55-i;
if  i>8&i<=29
midu=22994.15;
else
midu=7800;
end
%%%%%%%%%%%%%%%
其它部分的计算也解释一下，最好是上传word文档的公式几背景说明．
否则不是搞这个专业的高手也很难帮你修改程序．

yuanmm427

程序明显有问题！

hanxiao

当我运算完
yiii=q*i*p;
yrrr=q*j*p;
这两个命令之后就会出现下面的错误：

??? Error using ==> reshape
To RESHAPE the number of elements must not change.

Error in ==> sym.maple at 94
      result = reshape(result,size(varargin{3}));

Error in ==> sym.simplify at 14
r = maple('map','simplify',s);


另外还有：
（1）其实这个程序，当a=[z^2/68^2,z^3/68^3,z^4/68^4,z^5/68^5,z^6/68^6]时，其他的全不做修改，是完全可以运行正确的，而且结果也很正确，我和别的方法比较了。
（2）可是当a=[1-cos(0.0231*z),1-cos(0.0693*z),1-cos(0.115*z),1-cos(0.162*z),1-cos(0.208*z)]，即a 取三角函数的时候，就出现了这个错误。可能是化简除了问题，可是为什么化简有问题呢？

跪求各位大虾了！！

hanxiao

到底是怎么回事啊？郁闷阿。。。。

hanxiao

?

xjzuo

请先把问题和背景讲清楚一些.