为什么每次运行的图像都不一样呢？

vitihl

下面这段程序，为什么每次运行的图像都不一样呢？
clear all;
m1=3000; m2=2065; l=660; c1=5; k1=260; c2=5; k2=260; c3=30;  k3=150;
e1=0; e2=0; e3=189+e2; e4=189-e2; e11=l+e1; e12=l-e1; e21=l+e2; e22=l-e2;
w=44; j1=5.194e+6; j2=8.477e+6; F0=37.5; time=2*pi/w index=1;
M=zeros(1:1);%记录结果，是一个4维数组，分别记录了 A1 A2 A3 A4
for t=0:0.001:time,
   
a11=(-(e12*m1*w^2*sin(w*t)/2*l)+c1*w*cos(w*t)+k1*sin(w*t));
a12=(-(e11*m1*w^2*sin(w*t)/2*l)+c2*w*cos(w*t)+k2*sin(w*t));
a13=(-c1*w*cos(w*t)-k1*sin(w*t));
a14=(-c2*w*cos(w*t)-k2*sin(w*t));
a21=(-(j1*w^2*sin(w*t)/2*l)+c1*e11*w*cos(w*t)+k1*e11*sin(w*t));
a22=(-(j1*w^2*sin(w*t)/2*l)+c2*e12*w*cos(w*t)+k2*e12*sin(w*t));
a23=(-c1*e11*w*cos(w*t)-k1*e11*sin(w*t));
a24=(-c2*e12*w*cos(w*t)-k2*e12*sin(w*t));
a31=(-c1*w*cos(w*t)-k1*sin(w*t));
a32=(-c2*w*cos(w*t)-k2*sin(w*t));
a33=(-(e22*m2*w^2*sin(w*t)/2*l)+(c1+c3)*w*cos(w*t)+(k1+k3)*sin(w*t));
a34=(-(e21*m2*w^2*sin(w*t)/2*l)-(c2+c3)*w*cos(w*t)+(k2+k3)*sin(w*t));
a41=(-c1*e21*w*cos(w*t)-k1*e21*sin(w*t));
a42=(c2*e22*w*cos(w*t)+k2*e22*sin(w*t));
a43=(-(j2*w^2*sin(w*t)/2*l)+(c1+c3)*e21*w*cos(w*t)+(k1+k3)*e21*sin(w*t));
a44=((j2*w^2*sin(w*t)/2*l)-(c2+c3)*e22*w*cos(w*t)-(k2+k3)*e22*sin(w*t));
v=solve('a11*A1+a12*A2+a13*A3+a14*A4=0','a21*A1+a22*A2+a23*A3+a24*A4=0','a31*A1+a32*A2+a33*A3+a34*A4=2*F0*sin(w*t)','a41*A1+a42*A2+a43*A3+a44*A4=(e4-e3)*F0*sin(w*t)','A1,A2,A3,A4');
    %tsubs(v.A1),subs(v.A2),subs(v.A3),subs(v.A4)
    M(1,index)=subs(v.A1); M(2,index)=subs(v.A2);
    M(3,index)=subs(v.A3); M(4,index)=subs(v.A4);
    index=index+1;%自增长下标
end
plot(M(1,:),':')%画出A1的所有计算结果

[ 本帖最后由 ChaChing 于 2009-12-11 18:11 编辑 ]

sigma665

M=zeros(1:1);%

这个，是这样的吗

raozel

for t=0:0.001:time
建议改为for t=eps:0.001:time
0处求解易出奇异矩阵问题

vitihl

原帖由 raozel 于 2008-4-16 19:32 发表

                               
登录/注册后可看大图

for t=0:0.001:time
建议改为for t=eps:0.001:time
0处求解易出奇异矩阵问题

不行，还是每次运行结果都不一致

vitihl

对于之前在论坛里提到的那个四自由度微分方程组，编写如下程序
function dy=ep(t,y)
m1=3000; m2=2065; l=660; c1=5; k1=260; c2=5; k2=260; c3=30; k3=150;
e1=0; e2=0; e3=189+e2; e4=189-e2; e11=l+e1; e12=l-e1; e21=l+e2; e22=l-e2;
w=44; j1=530; j2=865; F0=37.5;
dy=[y(2);
    ((-(j1+e11^2*m1)*c1*y(2)-(j1-e11*e12*m1)*c2*y(4))+(j1+e11^2*m1)*c1*y(6)+(j1-e11*e12*m1)*c2*y(8)-(j1+e11^2*m1)*k1*y(1)-(j1-e11*e12*m1)*k2*y(3)+(j1+e11^2*m1)*k1*y(5)+(j1-e11*e12*m1)*k2*y(7))*2*l/((e11+e12)*m1*j1);
    y(4);
    ((-(j1+e11*e12*m1)*c1*y(2)-(j1-e12^2*m1)*c2*y(4))+(j1+e11*e12*m1)*c1*y(6)+(j1-e12^2*m1)*c2*y(8)-(j1+e11*e12*m1)*k1*y(1)-(j1-e12^2*m1)*k2*y(3)+(j1+e11*e12*m1)*k1*y(5)+(j1-e12^2*m1)*k2*y(7))*2*l/((e11-e12)*m1*j1);
    y(6);
    (j2+e21^2*m2)*c1*y(2)+(j2-e21*e22*m2)*c2*y(4)-(j2+e21^2*m2)*(c1+c3)*y(6)-(j2-e21*e22*m2)*(c2+c3)*y(8)+(j2+e21^2*m2)*k1*y(1)+(j2-e21*e22*m2)*k2*y(3)-(j2+e21^2*m2)*(k1+k3)*y(5)-(j2-e21*e22*m2)*(k2+k3)*y(7)+(2*j2+(e4-e3)*e21*m2)*F0*sin(w*t);
    y(8);
    (j2-e21*e22*m2)*c1*y(2)+(j2+e22^2*m2)*c2*y(4)-(j2-e21*e22*m2)*(c1+c3)*y(6)-(j2+e22^2*m2)*(c2+c3)*y(8)+(j2-e21*e22*m2)*k1*y(1)+(j2+e22^2*m2)*k2*y(3)-(j2-e21*e22*m2)*(k1+k3)*y(5)-(j2+e22^2*m2)*(k2+k3)*y(7)+(2*j2-(e4-e3)*e22*m2)*F0*sin(w*t)];

调用ode45
y0=[1 1 1 1 1 1 1 1]  
tic,[t,x]=ode45('ep',[0 0.2],y0);toc
  length(t),plot(x(:,1),x(:,3),x(:,5),x(:,7),':')
运行后x值显示NaN是否意味着方程无解

[ 本帖最后由 ChaChing 于 2009-12-11 18:17 编辑 ]

vitihl

Warning: Divide by zero.
(Type "warning off MATLAB:divideByZero" to suppress this warning.)
> In C:\MATLAB6p5\work\HouHouMatLab\ep.m at line 23
  In C:\MATLAB6p5\toolbox\matlab\funfun\private\odearguments.m at line 104
  In C:\MATLAB6p5\toolbox\matlab\funfun\ode45.m at line 155
  In C:\MATLAB6p5\work\HouHouMatLab\sy.m at line 2
程序报错内容

papayadong

老兄,我也是这个方程,解不出来!
你有什么办法?
帮下忙!
keyirou2000@eyou.com

[ 本帖最后由 eight 于 2008-4-22 10:53 编辑 ]

eight

原帖由 papayadong 于 2008-4-22 10:52 发表

                               
登录/注册后可看大图

老兄,我也是这个方程,解不出来!
你有什么办法?
帮下忙!
跪求!
keyirou2000@eyou.com

新手发帖请先看版规（存在不必要的字眼），警告一次

vitihl

function f=ep(t,y)
m1=3000; m2=2065; l=660; c1=5; k1=260; c2=5; k2=260; c3=30;  k3=150;
e1=0; e2=0; e3=189+e2; e4=189-e2; e11=l+e1; e12=l-e1; e21=l+e2; e22=l-e2;
w=44; j1=530; j2=865; F0=37.5;
f=[y(2);
    ((-(j1+e11^2*m1)*c1*y(2)-(j1-e11*e12*m1)*c2*y(4))+(j1+e11^2*m1)*c1*y(6)+(j1-e11*e12*m1)*c2*y(8)-(j1+e11^2*m1)*k1*y(1)-(j1-e11*e12*m1)*k2*y(3)+(j1+e11^2*m1)*k1*y(5)+(j1-e11*e12*m1)*k2*y(7))*2*l/((e11+e12)*m1*j1);
    y(4);
    ((-(j1+e11*e12*m1)*c1*y(2)-(j1-e12^2*m1)*c2*y(4))+(j1+e11*e12*m1)*c1*y(6)+(j1-e12^2*m1)*c2*y(8)-(j1+e11*e12*m1)*k1*y(1)-(j1-e12^2*m1)*k2*y(3)+(j1+e11*e12*m1)*k1*y(5)+(j1-e12^2*m1)*k2*y(7))*2*l/((e11-e12)*m1*j1);
    y(6);
    ((j2+e21^2*m2)*c1*y(2)+(j2-e21*e22*m2)*c2*y(4)-(j2+e21^2*m2)*(c1+c3)*y(6)-(j2-e21*e22*m2)*(c2+c3)*y(8)+(j2+e21^2*m2)*k1*y(1)+(j2-e21*e22*m2)*k2*y(3)-(j2+e21^2*m2)*(k1+k3)*y(5)-(j2-e21*e22*m2)*(k2+k3)*y(7)+(2*j2+(e4-e3)*e21*m2)*F0*sin(w*t))*2*l/((e21+e22)*m1*j1);
    y(8);
    ((j2-e21*e22*m2)*c1*y(2)+(j2+e22^2*m2)*c2*y(4)-(j2-e21*e22*m2)*(c1+c3)*y(6)-(j2+e22^2*m2)*(c2+c3)*y(8)+(j2-e21*e22*m2)*k1*y(1)+(j2+e22^2*m2)*k2*y(3)-(j2-e21*e22*m2)*(k1+k3)*y(5)-(j2+e22^2*m2)*(k2+k3)*y(7)+(2*j2-(e4-e3)*e22*m2)*F0*sin(w*t))*2*l/((e21+e22)*m1*j1)];

[ 本帖最后由 ChaChing 于 2009-12-11 18:18 编辑 ]