[求助]麻烦大家帮我看看这个程序～～

lujixx

clc;
clear all;
g=9.8;
c1=1.25*10^(-5);
delta=4.5852*10^(-2);
beta=1.1404*10^(-2);
omega=200;
omegag=omega*g^(1/2)/c1^(1/2);
t(1)=0;
h=0.1;
n=1;
C=[0 0 1 0;0 0 0 1;0 0 0 0;0 0 0 0];
x(:,1)=[0.0001;0.0004;0.0002;0.0006];
while t<=10
F=[0;0;beta*cos(omegag*t)-delta/omega*((2*x(1,n)^2+2*x(2,n)^2-4*x(1,n)*x(4,n)+...
4*x(2,n)*x(3,n))/((1-x(1,n)^2-x(2,n)^2)*(2+x(1,n)^2+x(2,n)^2))+(x(1,n)*x(3,n)+...
x(2,n)*x(4,n))*(pi-16/pi/(2+x(1,n)^2+x(2,n)^2))/(x(1,n)^2+x(2,n)^2)^(1/2)/(1-x(1,n)^2-x(2,n)^2)^(3/2))*cos(omegag*t)-...
delta/omega*(pi*(x(1,n)^2+x(2,n)^2)*(1-2*(x(1,n)*x(4,n)-x(2,n)*x(3,n))/(x(1,n)^2+x(2,n)^2))/(1-x(1,n)^2-x(2,n)^2)^(1/2)/(2+x(1,n)^2+...
x(2,n)^2)+4*(x(1,n)*x(3,n)+x(2,n)*x(4,n))/(2+x(1,n)^2+x(2,n)^2)/(1-x(1,n)^2-x(2,n)^2))*sin(omegag*t);
beta*sin(omegag*t)-1/omega^2-delta/omega*((2*x(1,n)^2+2*x(2,n)^2-4*x(1,n)*x(4,n)+4*x(2,n)*...
x(3,n))/((1-x(1,n)^2-x(2,n)^2)*(2+x(1,n)^2+x(2,n)^2))+(x(1,n)*x(3,n)+x(2,n)*x(4,n))*...
(pi-16/pi/(2+x(1,n)^2+x(2,n)^2))/(x(1,n)^2+x(2,n)^2)^(1/2)/(1-x(1,n)^2-x(2,n)^2)^(3/2))*...
sin(omegag*t)+delta/omega*(pi*(x(1,n)^2+x(2,n)^2)*(1-2*(x(1,n)*x(4,n)-x(2,n)*x(3,n))/..
(x(1,n)^2+x(2,n)^2))/(1-x(1,n)^2-x(2,n)^2)/(2+x(1,n)^2+x(2,n)^2)+4*(x(1,n)*x(3,n)+x(2,n)*x(4,n))/(2+x(1,n)^2+...
x(2,n)^2)/(1-x(1,n)^2-x(2,n)^2))*cos(omegag*t)];

k1=h*(C*x(:,n)+F);
k2=h*(C*(x(:,n)+k1/2)+F);
k3=h*(C*(x(:,n)+k2/2)+F);
k4=h*(C*(x(:,n)+k3)+F);
x(:,n+1)=x(:,n)+(1/6)*(k1+2*k2+2*k3+k4);
t(n+1)=n*h;
n=n+1;
end
x
plot(t,x(1,:));
我是个新手，以上是我自己弄的，也不知道问题在那儿。

happy

附件是空的

lujixx

不好意思～<BR>我已经改过了。<BR>象这种比较复杂的方程是不是不适合用ode45来算？

happy

<P>附件还是没内容</P>