[求助]常微分方程组求解，望高手帮看看到底哪里的问题，谢谢

WF1987

??? 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
[t,y]=ode45(@weifenzu,[0,0.2],x0,options);

WF1987

%M函数文件共五个，第一个
x0=[0.32;1500;3000;0.1;300;4;2500;293;0.000003];  %y（1）-y(9)的初值

options=odeset('relTol',1e-3,'Maxstep',0.00001);  %设置相对误差、最大积分步长

[t,y]=ode45(@weifenzu,[0,0.2],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('颗粒直径')

WF1987

%第二个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

WF1987

%第三个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  

WF1987

%第四个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

WF1987

%第五个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

WF1987

弄了两个多月了，都没出结果，导师催了，望高手帮忙看看，非常感谢

ChaChing

太多错误了!
1.函数与变数混用, 如f.m中Cd, 建议函数名多取些, 不要老是f, J,..., 容易出错不自知
2.函数参数个数不符, 如f.m有6个参数, 但出现仅3个参数呼叫

还有儘量简化程序, 式子那麼长容易出错不自知

WF1987

恩  好  我在重新写   回头再来

WF1987

我进行了修改，两个M文件如下：
x0=[0.32;1500;3000;0.1;300;4;2500;293;0.000003];  %y（1）-y(9)的初值

options=odeset('relTol',1e-3,'Maxstep',0.00001);  

[t,y]=ode45(@weifenzu,[0,0.2],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

%%可以运行，但没有正常的图形输出，什么情况呢？？ 谢谢

ChaChing

相信LZ可以更正相关错误, 毕竟自己程序自己最清楚!
若有结果记得分享, 当然再有问题, 也可再讨论, 虽说不见得能帮忙!:@)

WF1987

回复 11 # ChaChing 的帖子

谢谢  我再看看吧

WF1987

好久没来看看了