likefeng 发表于 2008-12-5 21:40

用ode45编程，程序有误，大家给解决一下。

用ode45编程，程序有误，大家给解决一下。
function dx=fz1(t,x,F1,F2,C12,C21,K11,K12,w1,w2)
dx=[x(2);
    0.002337*(F1-C12.*x(4)-K11*x(1)-K12*x(3))-0.02*w1*x(2)-w1^2*x(1);
    x(4);
0.00117357*(F2-C21*x(2)-K21*x(1)-K22*x(3))-0.02*w2*x(4)-w2^2*x(3)];





x0=;
syms M1 M2 A B C F1 F2 f1 f2 w t b w1;
syms s K1 K2 K3 K4K5 K11 K12 K21 K22;
s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;
T0=2*pi/w;t=0:T0/500:T0;
A=-0.3*w^2*cos(w*t).*cos(s*t)-s*(-0.3)*w*sin(s*t);
B=0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*cos(s*t);
C=-0.3*w^2*cos(s*t);
f1=48535.1+993.93*cos(s*t);
f2=-10426.4+(-177.19)*cos(s*t);
F1=f1+(-28510.61)*A.*b+(-6619.84)*B.*b^2+(-872.82)*C;
F2=f2+468.19*A*b+(-47391.25)*B*b^2+75.82*C;
C12=1078.62*(-0.3)*b*sin(s*t);
C21=1078.62*0.3*b*sin(s*t);
K1=b^2*0.09*w^(2)*sin(w*t).*sin(w*t);
K2=b*((-0.3)*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K3=b*(0.3*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K4=-b*(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*sin(s*t)-s^2;
K5=b^(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*cos(s*t).*cos(s*t);
K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;
K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;
K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;
K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;
=ode45('fz1',t,x0);
subplot(2,2,1);
plot(t,x(:,1));
hold on;
subplot(2,2,2);
plot(t,x(:,2));
hold on;
subplot(2,2,3);
plot(t,x(:,3));
hold on;
subplot(2,2,4);
plot(t,x(:,4));


Error in ==> fz1 at 2
dx=[x(2);
Error in ==> funfun\private\odearguments at 110
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> zfz1 at 26
=ode45('fz1',t,x0);

likefeng 发表于 2008-12-6 10:37

问题已解决，正确如下：
function dx=gfz1(t,x,F1,F2,C12,C21,K11,K12,K21,K22,w1,w2)
s=1.6;b=1;w=4.8;w1=5.14;w2=9.49;
A=-0.3*w^2*cos(w*t).*cos(s*t)-s*(-0.3)*w*sin(s*t);
B=0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*cos(s*t);
C=-0.3*w^2*cos(s*t);
f1=48535.1+993.93*cos(s*t);
f2=-10426.4+(-177.19)*cos(s*t);
F1=f1+(-28510.61)*A.*b+(-6619.84)*B.*b^2+(-872.82)*C;
F2=f2+468.19*A*b+(-47391.25)*B*b^2+75.82*C;
C12=1078.62*(-0.3)*b*sin(s*t);
C21=1078.62*0.3*b*sin(s*t);
K1=b^2*0.09*w^(2)*sin(w*t).*sin(w*t);
K2=b*((-0.3)*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K3=b*(0.3*w^(2)*cos(w*t).*sin(s*t)-0.3*w*sin(w*t).*cos(s*t)*s);
K4=-b*(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*sin(s*t).*sin(s*t)-s^2;
K5=b^(2)*0.09*w^(2)*sin(w*t).*sin(w*t).*cos(s*t).*cos(s*t);
K11=381.34*K1+81.22*K2+81.22*K3+17.72*K4+806.29*K5;
K12=-115.92*K1+517.92*K2+(-21.4)*K3+115.95*K4+(-206.29)*K5;
K21=-115.92*K1-21.42*K2+517.92*K3+115.95*K4+(-206.29)*K5;
K22=61.04*K1-114.93*K2-114.93*K3+784.94*K4+1373.31*K5;
dx=[x(2);
    0.002337*(F1-C12*x(4)-K11*x(1)-K12*x(3))-0.02*w1*x(2)-w1^(2)*x(1);
    x(4);
0.00117357*(F2-C21*x(2)-K21*x(1)-K22*x(3))-0.02*w2*x(4)-w2^(2)*x(3)];

T0=2*pi/w;t=0:T0/500:T0;
x0=;
=ode45('gfz1',t,x0);

likefeng 发表于 2008-12-6 10:39

另外想请教大家，根据上面的微分方程如何得到附件图，图中W0就是程序中的W。

gwcrow 发表于 2008-12-10 11:25

求助：求解微分方程

dx/dt=x(1-Cx)-Nx(1+x),c<1,N在0.5至1.5之间，用ode45,请问怎么编程，本人刚刚接触matlab，现在弄得一头雾水，还有请教ode45使用的格式，万分感谢！！

likefeng 发表于 2008-12-10 15:11

找本matlab的书都有例题，或论坛。
页: [1]
查看完整版本: 用ode45编程，程序有误，大家给解决一下。