声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1841|回复: 5

[编程技巧] 请教程序中reshape出错

[复制链接]
发表于 2007-4-25 18:23 | 显示全部楼层 |阅读模式

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

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

x
??? 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.int at 51
   r = reshape(maple('map','int',f(:),[x.s '=(' a.s ')..(' b.s ')']),size(f));

Error in ==> fgxczdlxwt at 39
    k2=int(ee*ii*(diff(ph2,2))'*(diff(ph2,2)),0,l);
我的原程序如下:
%用newmark法求解非惯性场中的振动问题
clc
clear
%定义柔性梁材料参数,l为柔性梁的长度,a为横截面积,ii为转动惯量
l=0.9;  aa=3.18*1e-5;   ii=2.65*1e-12;
p=7866; ee=2.01*1e11;  r0=0.55;
n=input('请输入所取的模态阶数=')    %%%  n是取的阶数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=0.05;
t_0=0;
t_max=100;
jisuanbushu=(t_max-t_0)/dt+1;
u=zeros(2*n,jisuanbushu);
v=zeros(2*n,jisuanbushu);
a=zeros(2*n,jisuanbushu);
%计算模态矩阵ph1为纵向变形的,ph2为横向变形的         
syms x t y
for i=1:n
    ph1(i)=sin((2*i-1)*pi*x/(2*l));
    if i==1
        ph2(i)=cos(1.875*x/l)-cosh(1.875*x/l)-(cos(1.875)+cosh(1.875))*(sin(1.875*x/l)-sinh(1.875*x/l))/(sin(1.875)+sinh(1.875));
    elseif i==2
        ph2(i)=cos(4.694*x/l)-cosh(4.694*x/l)-(cos(4.694)+cosh(4.694))*(sin(4.694*x/l)-sinh(4.694*x/l))/(sin(4.694)+sinh(4.694));
    else
        ph2(i)=cos((i-0.5)*pi*x/l)-cosh((i-0.5)*pi*x/l)-(cos((i-0.5)*pi)+cosh((i-0.5)*pi))*(sin((i-0.5)*pi*x/l)-sinh((i-0.5)*pi*x/l))/(sin((i-0.5)*pi)+sinh((i-0.5)*pi));
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=0:1:jisuanbushu
    t=0;
    if  t<=80
        omega=4*(1-exp(-5*t/60));
    else
        omega=4;
    end
    %计算各系数矩阵的值
    clear maplemex
    k1=int(ee*aa*(diff(ph1))'*(diff(ph1)),0,l);
    k2=int(ee*ii*(diff(ph2,2))'*(diff(ph2,2)),0,l);
    u01=int(p*aa*ph1,0,l);
    u02=int(p*aa*ph2,0,l);
    u11=int(p*aa*x*ph1,0,l);
    u12=int(p*aa*x*ph2,0,l);
    s=int((diff(subs(ph2,x,y)))'*(diff(subs(ph2,x,y))),y,0,x);
    d0=int(p*aa*s,0,l);
    d1=int(p*aa*x*s,0,l);
    rr=int(p*aa*ph1'*ph2,0,l);

    m1=int(p*aa*ph1'*ph1,0,l);
    m2=int(p*aa*ph2'*ph2,0,l);
    mq1c=-rr*q2;
    mq2c=(r0*u02+u12+q1'*rr)';
    gq1q2=-rr;
    gq2q1=rr';
    kq1q1=k1-(omega)^2*m1;
    kq2q2=k2-(omega)*m2+(omega)^2*(r0*d0+d1);
    qq1=(omega)^2*(r0*u01'+u11');

    m=[m1,0;0,m2];
    c=2*omega*[0,gq1q2;gq2q1,0];
    k=[kq1q1,0;0,kq2q2];
    q=[qq1-(diff(omega))*mq1c;mq2c];


    %利用newmark法求解梁的振动位移
    b=0.25;
    r=0.5;

    a0=1/(b*(dt)^2);
    a1=r/(b*dt);
    a2=1/(b*dt);
    a3=1/(2*b)-1;
    a4=r/b-1;
    a5=0.5*dt*((r/b)-2);
    a6=dt*(1-r);
    a7=r*dt;


    k=k+a0*m+a1*c;
    q=q+m*(a6*u(:,i)+a2*v(:,i)+a3*a(:,i))+c*(a1*u(:,i)+a4*v(:,i)+a5*a(:,i));
    u(:,i+1)=inv(k)*q(:,i+1);
    v(:,i+1)=a0*(x(:,i+1)-x(:,i))-a2*dx(:,i)-a3*ddx(:,i);
    a(:,i+1)=v(:,i)+a6*a(:,i)+a7*a(:,i+1);
end
回复
分享到:

使用道具 举报

 楼主| 发表于 2007-4-25 18:26 | 显示全部楼层

回复 #1 satlxl 的帖子

我看到有人关于这类帖子,
Y=vpa(Y,4);   %其中4可以根据需要调整;
    if mod(n,500)==0   %其中500可以根据需要调整,越大越好
       clear maplemex
    end
但是没看懂,请哪位大虾帮我修改一下,非常感谢
发表于 2007-4-25 18:32 | 显示全部楼层
原帖由 satlxl 于 2007-4-25 18:23 发表
??? 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.int at 51
  ...



reshape 函数要求变换前后的矩阵总大小不变
 楼主| 发表于 2007-4-25 18:38 | 显示全部楼层
这个我知道,但是怎么改呢?我只是对一个1*n的矩阵求导转置,然后积分,应该没改变大小阿
发表于 2007-4-25 18:40 | 显示全部楼层
原帖由 satlxl 于 2007-4-25 18:26 发表
我看到有人关于这类帖子,
Y=vpa(Y,4);   %其中4可以根据需要调整;
    if mod(n,500)==0   %其中500可以根据需要调整,越大越好
       clear maplemex
    end
但是没看懂,请哪位大虾帮我修改一下,非 ...


看帖子:【已解决】计算过程中整数太大
 楼主| 发表于 2007-4-26 14:34 | 显示全部楼层
还是看不懂啊,不知道有没有哪本书介绍的详细一点阿,clear maplemex我也用了,还是不行,能就我的程序给点指点马?太感谢你了,急望回复!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-28 17:40 , Processed in 0.077498 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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