[求助]常微分方程组求解,望高手帮看看到底哪里的问题,谢谢
??? Attempted to access Cd(9); index out of bounds because numel(Cd)=1.Error in ==> f at 5
f=(3*x1*x4*Cd(Re)*(x5-x2)*abs(x5-x2))/(1.2*x9);
Error in ==> weifenzu at 14
dx=[(-f(x(1),x(4),Cd(Re),x(5),x(2),x(9))-J(x(8),x(1),Re)*x(2))/x(2)^2
Error in ==> 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 ==> ode15s at 5
=ode45(@weifenzu,,x0,options);
%M函数文件共五个,第一个
x0=;%y(1)-y(9)的初值
options=odeset('relTol',1e-3,'Maxstep',0.00001);%设置相对误差、最大积分步长
=ode45(@weifenzu,,x0,options); %调用M文件weifenzu。
subplot(3,3,1);plot(t,y(:,1),'k') %将当前窗口分割成3行3列区域,指定第1个编号区域为当前绘图区域。
xlabel('x(m)');ylabel('σp(kg/m^3)');title('颗粒浓度分布')%x、y坐标的名称及图形名称
subplot(3,3,2);plot(t,y(:,2),'r':t:y(:,5),'r')
xlabel('x(m)');ylabel('v(m/s)');title('速度分布')
subplot(3,3,3);plot(t,y(:,3),'g')
xlabel('x(m)');ylabel('');title('颗粒内能')
subplot(3,3,4);plot(t,y(:,4),'k')
xlabel('x(m)');ylabel('');title('气体浓度')
subplot(3,3,5);plot(t,y(:,6),'g')
xlabel('x(m)');ylabel('P(MPa)');title('气体压力')
subplot(3,3,6);plot(t,y(:,7),'k-')
xlabel('x(m)');ylabel('T(K)');title('气体温度')
subplot(3,3,7);plot(t,y(:,9),'k-')
xlabel('x(m)');ylabel('d(um)');title('颗粒直径') %第二个M文件
function dx=weifenzu(t,x)
u=1; P=1;L=1;E=1;B=1;D=1652;F=1; %先假设各参量都为常数1,后期再进行修正
Re=x(4)*abs(x(5)-x(2))*x(9)/u; %雷诺数的表达式
q=(12+2.754*(x(4)*abs(x(5)-x(2))*x(9)/u)^0.55*P^0.33)*(x(7)-x(8))*L*x(1)/x(9)^2+6*x(1)*(x(7)^4-x(8)^4)*E*B/x(9);
dx=[(-f(x(1),x(4),Cd(Re),x(5),x(2),x(9))-J(x(8),x(1),Re)*x(2))/x(2)^2
f(x(1),x(4),Cd(Re),x(5),x(2),x(9))/(x(1)*x(2))
q/(x(1)*x(2))
((-0.4*D^3+2.2*x(5)*D^2+36.256*x(7)*D-10623*D-3.2*D*x(5)^2+90.64*x(5)*x(7)+10496.112*x(5)+1.4*x(5)^3)*0.046*x(4)^2+(-2.4*J(x(8),x(1),Re)*x(5)^2+1.4*J(x(8),x(1),Re)*x(2)*x(5)-0.4*J(x(8),x(1),Re)*x(3)-1.4*f(x(1),x(4),Cd(Re))*x(5)+0.4*f(x(1),x(4),Cd(Re))*x(2)+0.4*q)*x(4)+1.4*J(x(8),x(1),Re)*x(6))/(1.4*x(5)*x(6)-x(4)*x(5)^3)
(-0.018*x(4)*D^3+0.101*x(4)*x(5)*D^2-488.658*x(4)*D+1.668*x(4)*x(7)*D-0.147*x(4)*D*x(5)^2+488.658*x(4)*x(5)-1.668*x(4)*x(5)*x(7)+0.064*x(4)*x(5)^3+0.4*q-1.4*J(x(8),x(1),Re)*x(5)^2-0.4*J(x(8),x(1),Re)*x(3)+1.4*J(x(8),x(1),Re)*x(2)*x(5)+0.4*f(x(1),x(4),Cd(Re))*x(2)-1.4*f(x(1),x(4),Cd(Re))*x(5))/(x(4)*x(5)^2-1.4*x(6))
(0.018*x(5)^4-0.055*D*x(5)^3-(1.668*x(7)-0.055*D^2-488.658)*x(5)^2-488.658*D*x(5)-0.018*x(5)*D^3+1.668*x(5)*x(7)*D)*x(4)^2+(-0.4*J(x(8),x(1),Re)*x(5)^3+0.064*x(6)*x(5)^2+0.4*J(x(8),x(1),Re)*x(2)*x(5)^2-0.4*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(5)^2-0.129*D*x(5)*x(6)-0.4*J(x(8),x(1),Re)*x(3)*x(5)+0.4*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(2)*x(5)+0.4*q*x(5)+0.064*x(6)*D^2)*x(4)-1.4*J(x(8),x(1),Re)*x(5)*x(6)-1.4*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(6)+1.4*J(x(8),x(1),Re)*x(2)*x(6)
(0.002*x(4)^2*x(5)^5-(0.007*D*x(4)^2+0.048*J(x(8),x(1),Re)*x(4))*x(5)^4+(58.775*x(4)^2-0.265*x(7)*x(4)^2+0.007*D^2*x(4)^2-0.048*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(4)+0.048*J(x(8),x(1),Re)*x(2)*x(4)+0.008*x(6)*x(4))*x(5)^3+((-205.708*D+0.7*x(7)*D+0.147*x(7)*D-0.002*D^3+146.924*D-0.5*D*x(7))*x(4)^2+0.048*q*x(4)-0.048*J(x(8),x(1),Re)*x(3)*x(4)-0.015*D*x(6)*x(4)+2.4*J(x(8),x(1),Re)*x(4)*x(7)+0.048*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(2)*x(4)-0.168*J(x(8),x(1),Re)*x(6))*x(5)^2+(-488.658*x(7)*x(4)^2+1.668*x(4)^2*x(7)^2-0.101*x(7)*x(4)^2*D^2+0.008*x(6)*x(4)*D^2+1.4*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(7)*x(4)-1.4*J(x(8),x(1),Re)*x(2)*x(4)*x(7)+0.168*J(x(8),x(1),Re)*x(2)*x(6)-0.168*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(6))*x(5)+(1710.304*x(7)*D-5.837*D*x(7)^2+0.018*x(7)*D^3-1221.646*x(7)*D+4.169*D*x(7)^2)*x(4)^2+(0.4*J(x(8),x(1),Re)*x(3)*x(7)-0.4*q*x(7)-0.4*f(x(1),x(4),Cd(Re),x(5),x(2),x(9))*x(2)*x(7))*x(4)-1.4*J(x(8),x(1),Re)*x(6)*x(7))/(1.4*x(4)*x(5)*x(6)-x(4)^2*x(5)^3)
((x(7)-x(8))*(x(1)+x(4))^0.5)/(0.014*F^0.9*x(6)^0.5);
47.619*F^0.9+13.143*F^0.9*(x(4)*abs(x(5)-x(2))*x(9)/u)^0.5
];
end %第三个M文件
function Cd=Cd(Re)
if Re<800
Cd=24*(1+0.15*Re^0.687)/Re;
else Cd=24*(1+0.15*Re^0.687)/Re+0.42/(1+42500*Re^(-1.16));
end
end %第四个M文件
function f=f(x1,x4,Cd,x5,x2,x9)
u=1;
Re=25000*x4*abs(x5-x2)*x9/u;
f=(3*x1*x4*Cd(Re)*(x5-x2)*abs(x5-x2))/(1.2*x9);
end %第五个M文件
function J=J(x8,x1,Re)
if x8<933
J=0;
else J=142.857*x1*(1+0.276*sqrt(Re))*F^0.9;
end
end 弄了两个多月了,都没出结果,导师催了,望高手帮忙看看,非常感谢 太多错误了!
1.函数与变数混用, 如f.m中Cd, 建议函数名多取些, 不要老是f, J,..., 容易出错不自知
2.函数参数个数不符, 如f.m有6个参数, 但出现仅3个参数呼叫
还有儘量简化程序, 式子那麼长容易出错不自知
恩好我在重新写 回头再来 我进行了修改,两个M文件如下:
x0=;%y(1)-y(9)的初值
options=odeset('relTol',1e-3,'Maxstep',0.00001);
=ode45(@weifenzu,,x0,options); %调用M文件weifenzu。
subplot(3,3,1);plot(t,y(:,1),'k') %将当前窗口分割成3行3列区域,指定第1个编号区域为当前绘图区域。
xlabel('x(m)');ylabel('σp(kg/m^3)');title('颗粒浓度分布')%x、y坐标的名称及图形名称
subplot(3,3,2);plot(t,y(:,2),'r',t,y(:,5),'r')
xlabel('x(m)');ylabel('v(m/s)');title('速度分布')
subplot(3,3,3);plot(t,y(:,3),'g')
xlabel('x(m)');ylabel('');title('颗粒内能')
subplot(3,3,4);plot(t,y(:,4),'k')
xlabel('x(m)');ylabel('');title('气体浓度')
subplot(3,3,5);plot(t,y(:,6),'g')
xlabel('x(m)');ylabel('P(MPa)');title('气体压力')
subplot(3,3,6);plot(t,y(:,7),'k-')
xlabel('x(m)');ylabel('T(K)');title('气体温度')
subplot(3,3,7);plot(t,y(:,9),'k-')
xlabel('x(m)');ylabel('d(um)');title('颗粒直径')
%…………………………………………………………………………另外一个M文件如下:
function dx=weifenzu(t,x)
u=1; P=1;L=1;E=1;B=1;D=1652;F=1; %先假设各参量都为常数1,后期再进行修正
Re=x(4)*abs(x(5)-x(2))*x(9)/u; %雷诺数的表达式
q=(12+2.754*(x(4)*abs(x(5)-x(2))*x(9)/u)^0.55*P^0.33)*(x(7)-x(8))*L*x(1)/x(9)^2+6*x(1)*(x(7)^4-x(8)^4)*E*B/x(9);
CdRe=24*(1+0.15*Re^0.687)/Re*(Re<800)+(Re>=800)*(24*(1+0.15*Re^0.687)/Re+0.42/(1+42500*Re^(-1.16)));%0.42/(1+42500*Re^(-1.16)))
fF=(3*x(1)*x(4)*CdRe*(x(5)-x(2))*abs(x(5)-x(2)))/(1.2*x(9));
J=(x(8)>=933)*142.857*x(1)*(1+0.276*sqrt(Re))*F^0.9;
dx=[(-fF-J*x(2))/x(2)^2
fF/(x(1)*x(2))%fF/(x(1)*x(2))
q/(x(1)*x(2))
((-0.4*D^3+2.2*x(5)*D^2+36.256*x(7)*D-10623*D-3.2*D*x(5)^2+90.64*x(5)*x(7)+10496.112*x(5)+1.4*x(5)^3)*0.046*x(4)^2+(-2.4*J*x(5)^2+1.4*J*x(2)*x(5)-0.4*J*x(3)-1.4*fF*x(5)+0.4*fF*x(2)+0.4*q)*x(4)+1.4*J*x(6))/(1.4*x(5)*x(6)-x(4)*x(5)^3)
(-0.018*x(4)*D^3+0.101*x(4)*x(5)*D^2-488.658*x(4)*D+1.668*x(4)*x(7)*D-0.147*x(4)*D*x(5)^2+488.658*x(4)*x(5)-1.668*x(4)*x(5)*x(7)+0.064*x(4)*x(5)^3+0.4*q-1.4*J*x(5)^2-0.4*J*x(3)+1.4*J*x(2)*x(5)+0.4*fF*x(2)-1.4*fF*x(5))/(x(4)*x(5)^2-1.4*x(6))
(0.018*x(5)^4-0.055*D*x(5)^3-(1.668*x(7)-0.055*D^2-488.658)*x(5)^2-488.658*D*x(5)-0.018*x(5)*D^3+1.668*x(5)*x(7)*D)*x(4)^2+(-0.4*J*x(5)^3+0.064*x(6)*x(5)^2+0.4*J*x(2)*x(5)^2-0.4*fF*x(5)^2-0.129*D*x(5)*x(6)-0.4*J*x(3)*x(5)+0.4*fF*x(2)*x(5)+0.4*q*x(5)+0.064*x(6)*D^2)*x(4)-1.4*J*x(5)*x(6)-1.4*fF*x(6)+1.4*J*x(2)*x(6)
(0.002*x(4)^2*x(5)^5-(0.007*D*x(4)^2+0.048*J*x(4))*x(5)^4+(58.775*x(4)^2-0.265*x(7)*x(4)^2+0.007*D^2*x(4)^2-0.048*fF*x(4)+0.048*J*x(2)*x(4)+0.008*x(6)*x(4))*x(5)^3+((-205.708*D+0.7*x(7)*D+0.147*x(7)*D-0.002*D^3+146.924*D-0.5*D*x(7))*x(4)^2+0.048*q*x(4)-0.048*J*x(3)*x(4)-0.015*D*x(6)*x(4)+2.4*J*x(4)*x(7)+0.048*fF*x(2)*x(4)-0.168*J*x(6))*x(5)^2+(-488.658*x(7)*x(4)^2+1.668*x(4)^2*x(7)^2-0.101*x(7)*x(4)^2*D^2+0.008*x(6)*x(4)*D^2+1.4*fF*x(7)*x(4)-1.4*J*x(2)*x(4)*x(7)+0.168*J*x(2)*x(6)-0.168*fF*x(6))*x(5)+(1710.304*x(7)*D-5.837*D*x(7)^2+0.018*x(7)*D^3-1221.646*x(7)*D+4.169*D*x(7)^2)*x(4)^2+(0.4*J*x(3)*x(7)-0.4*q*x(7)-0.4*fF*x(2)*x(7))*x(4)-1.4*J*x(6)*x(7))/(1.4*x(4)*x(5)*x(6)-x(4)^2*x(5)^3)
((x(7)-x(8))*(x(1)+x(4))^0.5)/(0.014*F^0.9*x(6)^0.5);
47.619*F^0.9+13.143*F^0.9*(x(4)*abs(x(5)-x(2))*x(9)/u)^0.5
];
end
%%可以运行,但没有正常的图形输出,什么情况呢?? 谢谢 相信LZ可以更正相关错误, 毕竟自己程序自己最清楚!
若有结果记得分享, 当然再有问题, 也可再讨论, 虽说不见得能帮忙!:@) 回复 11 # ChaChing 的帖子
谢谢我再看看吧 好久没来看看了
页:
[1]