请教：计算结果的化简问题？

hanxiao

本人初学 matlab，遇到一个问题，困扰了我 3天了 ，怎么试也不行 ，请各位大侠帮帮忙阿！谢谢了！

问题是这样的：
计算结果是个很大的式子，我想把它化简，用simplify这个命令总是出错误，只能用vpa这个命令，结果得到的结果还是很大。请问，simplify和vpa这两个命令具体有什么区别阿？能不能举例说明。
谢谢了先！

pengweicai

函数  simplify
格式  R = simplify(S)
说明  使用Maple软件中的化简规则，将化简符号矩阵S中每一元素。

函数  vpa
格式  R = vpa(A)   %用可变精度算法来计算A中的每一元素，使其成为有d位精确度的十进制数。其中d为命令digits设置的当前位数。R中的每一元素为一符号表达式。
R = vpa(A,d)或R = vpa A d   %用参量d指定的位数(而非命令digits设置的位数)来表示A中的每一元素。R中的每一元素为一符号表达式。


  最好把错误 的提示 ，贴出来。

hanxiao

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

整个这个问题已经困扰我足足一个月了,各位大侠帮帮忙吧!谢谢了先!
以下是我的程序,帮我看看,万分感谢!


%-- 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)

xjzuo

请将背景知识稍微描述一下，把问题先讲清楚，否则很容易进行错误计算．
例如，从你的程序来看,似乎下面语句可以提出第一个循环,仅取i=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);
%%%%%%%%%
否则，应当加下标保存b0,b1,...

hanxiao

原帖由 hanxiao 于 2006-12-4 15:55 发表
(1)程序不复杂,计算结果十分复杂.
s38uu,s3801和s3802是我想要的计算结果,长得吓人!我想怎么能化简一下?
我是临时抱佛脚,刚刚学的matlab.想问一下 :
(2)对s38uu怎么化出它 的图形,并显示图形的数据结果?
(3) ...

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

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

hanxiao

跪求各位大虾帮帮忙啊