iamwho_2001
2010-09-24, 05:35
我想用下面的函数求解ODE， 但是总出现这个问题：
??? Error using ==> funfun\private\odearguments at 117
Solving VDP1 requires an initial condition vector of length 24.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> Untitled3 at 1
[t,y] = ode45('vdp1',[0 24],[10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0]);
但是我总也解决不了，所以请求大家的帮忙，多谢！
下面是函数和主程序：
子函数：
function out1 = vdp1(t,y)
%out1 = [y(2); (1-y(1)^2)*y(2) - y(1)];
km_minus=0.0002;
phT=5000;
kx1=2;
kmlat=5000;
km_plus=0.0002;
Ind3=370;
Ind4=1111;
Ind5=3333;
Ind6=10000;
Ind7=3333;
Ind8=1111;
kn_plus=130;
kx3=1300;
kmmpk=1000;
kx2=0.000002;
kn_minus=0.016;
out1 = [(km_minus*phT*y(2)) +(kx1*y(3)^2*y(2)/(kmlat^2+y(3)^2)) -(km_plus*y(1)*Ind3);
(km_plus*y(1)*Ind3) -(km_minus*phT*y(2)) -(kx1*y(3)^2*y(2)/(kmlat^2+y(3)^2));
(kn_plus) +(kx3*y(4)/(kmmpk*kmmpk+y(4))) -(kx2*y(3)*y(2)) -(kn_minus*y(3));
(y(6)*0.5)^2;
(kx1*y(7)^2*y(6)/(kmlat*kmlat+y(7)^2)) +(km_minus*phT*y(6)) -(km_plus*y(5)*Ind4);
-(kx1*y(7)^2*y(6)/(kmlat*kmlat+y(7)^2)) -(km_minus*phT*y(6)) +(km_plus*y(5)*Ind4);
(kx3*y(8)/(kmmpk*kmmpk+y(8))) +(kn_plus) -(kn_minus*y(7)) -(kx2*y(7)*y(6));
(y(2)+y(10)*0.5)^2;
(kx1*y(11)^2*y(10)/(kmlat*kmlat+y(11)^2)) +(km_minus*phT*y(10)) -(km_plus*y(9)*Ind4);
-(kx1*y(11)^2*y(10)/(kmlat*kmlat+y(11)^2)) -(km_minus*phT*y(10)) +(km_plus*y(9)*Ind4);
(kx3*y(12)/(kmmpk*kmmpk+y(12))) +(kn_plus) -(kn_minus*y(11)) -(kx2*y(11)*y(10));
(y(6)*0.5+y(14)*0.5)^2;
(kx1*y(15)^2*y(14)/(kmlat*kmlat+y(15)^2)) +(km_minus*phT*y(14)) -(km_plus*y(13)*Ind4);
-(kx1*y(15)^2*y(14)/(kmlat*kmlat+y(15)^2)) -(km_minus*phT*y(14)) +(km_plus*y(13)*Ind4);
(kx3*y(16)/(kmmpk*kmmpk+y(16))) +(kn_plus) -(kn_minus*y(15)) -(kx2*y(15)*y(14));
(y(10)*0.5+y(18)*0.5)^2;
(kx1*y(19)^2*y(18)/(kmlat*kmlat+y(19)^2)) +(km_minus*phT*y(18)) -(km_plus*y(17)*Ind4);
-(kx1*y(19)^2*y(18)/(kmlat*kmlat+y(19)^2)) -(km_minus*phT*y(18)) +(km_plus*y(17)*Ind4);
(kx3*y(20)/(kmmpk*kmmpk+y(20))) +(kn_plus) -(kn_minus*y(19)) -(kx2*y(19)*y(18));
(y(14)*0.5+y(22))^2;
(kx1*y(23)^2*y(22)/(kmlat*kmlat+y(23)^2)) +(km_minus*phT*y(22)) -(km_plus*y(21)*Ind4);
-(kx1*y(23)^2*y(22)/(kmlat*kmlat+y(23)^2)) -(km_minus*phT*y(22)) +(km_plus*y(21)*Ind4);
(kx3*y(24)/(kmmpk*kmmpk+y(24))) +(kn_plus) -(kn_minus*y(23)) -(kx2*y(23)*y(22));
(y(18)*0.5)^2]
end
主程序：
[t,y] = ode45('vdp1',[0 24],[10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0]);
%plot(t,y(:,1),'-',t,y(:,2),'-.')
??? Error using ==> funfun\private\odearguments at 117
Solving VDP1 requires an initial condition vector of length 24.
Error in ==> ode45 at 173
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in ==> Untitled3 at 1
[t,y] = ode45('vdp1',[0 24],[10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0]);
但是我总也解决不了，所以请求大家的帮忙，多谢！
下面是函数和主程序：
子函数：
function out1 = vdp1(t,y)
%out1 = [y(2); (1-y(1)^2)*y(2) - y(1)];
km_minus=0.0002;
phT=5000;
kx1=2;
kmlat=5000;
km_plus=0.0002;
Ind3=370;
Ind4=1111;
Ind5=3333;
Ind6=10000;
Ind7=3333;
Ind8=1111;
kn_plus=130;
kx3=1300;
kmmpk=1000;
kx2=0.000002;
kn_minus=0.016;
out1 = [(km_minus*phT*y(2)) +(kx1*y(3)^2*y(2)/(kmlat^2+y(3)^2)) -(km_plus*y(1)*Ind3);
(km_plus*y(1)*Ind3) -(km_minus*phT*y(2)) -(kx1*y(3)^2*y(2)/(kmlat^2+y(3)^2));
(kn_plus) +(kx3*y(4)/(kmmpk*kmmpk+y(4))) -(kx2*y(3)*y(2)) -(kn_minus*y(3));
(y(6)*0.5)^2;
(kx1*y(7)^2*y(6)/(kmlat*kmlat+y(7)^2)) +(km_minus*phT*y(6)) -(km_plus*y(5)*Ind4);
-(kx1*y(7)^2*y(6)/(kmlat*kmlat+y(7)^2)) -(km_minus*phT*y(6)) +(km_plus*y(5)*Ind4);
(kx3*y(8)/(kmmpk*kmmpk+y(8))) +(kn_plus) -(kn_minus*y(7)) -(kx2*y(7)*y(6));
(y(2)+y(10)*0.5)^2;
(kx1*y(11)^2*y(10)/(kmlat*kmlat+y(11)^2)) +(km_minus*phT*y(10)) -(km_plus*y(9)*Ind4);
-(kx1*y(11)^2*y(10)/(kmlat*kmlat+y(11)^2)) -(km_minus*phT*y(10)) +(km_plus*y(9)*Ind4);
(kx3*y(12)/(kmmpk*kmmpk+y(12))) +(kn_plus) -(kn_minus*y(11)) -(kx2*y(11)*y(10));
(y(6)*0.5+y(14)*0.5)^2;
(kx1*y(15)^2*y(14)/(kmlat*kmlat+y(15)^2)) +(km_minus*phT*y(14)) -(km_plus*y(13)*Ind4);
-(kx1*y(15)^2*y(14)/(kmlat*kmlat+y(15)^2)) -(km_minus*phT*y(14)) +(km_plus*y(13)*Ind4);
(kx3*y(16)/(kmmpk*kmmpk+y(16))) +(kn_plus) -(kn_minus*y(15)) -(kx2*y(15)*y(14));
(y(10)*0.5+y(18)*0.5)^2;
(kx1*y(19)^2*y(18)/(kmlat*kmlat+y(19)^2)) +(km_minus*phT*y(18)) -(km_plus*y(17)*Ind4);
-(kx1*y(19)^2*y(18)/(kmlat*kmlat+y(19)^2)) -(km_minus*phT*y(18)) +(km_plus*y(17)*Ind4);
(kx3*y(20)/(kmmpk*kmmpk+y(20))) +(kn_plus) -(kn_minus*y(19)) -(kx2*y(19)*y(18));
(y(14)*0.5+y(22))^2;
(kx1*y(23)^2*y(22)/(kmlat*kmlat+y(23)^2)) +(km_minus*phT*y(22)) -(km_plus*y(21)*Ind4);
-(kx1*y(23)^2*y(22)/(kmlat*kmlat+y(23)^2)) -(km_minus*phT*y(22)) +(km_plus*y(21)*Ind4);
(kx3*y(24)/(kmmpk*kmmpk+y(24))) +(kn_plus) -(kn_minus*y(23)) -(kx2*y(23)*y(22));
(y(18)*0.5)^2]
end
主程序：
[t,y] = ode45('vdp1',[0 24],[10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0; 10000; 0; 0; 0]);
%plot(t,y(:,1),'-',t,y(:,2),'-.')