24小时热门版块排行榜

返回列表

小木虫zb

木虫 (正式写手)

应助: 11 (小学生)
金币: 1560.9
红花: 5
帖子: 345
在线: 142.1小时
虫号: 1227944
注册: 2011-03-10
性别: GG
专业: 分离过程

[求助] matlab程序的问题

小弟用matlab计算一个偏微分方程，程序运行不出来，求大神相助，程序如下
function DBT_MNA%DBT用1表示,MNA用2表示
clear all;clc
A=importdata('DBT+MNA-RS-03.txt');%数据源文件中横坐标为g-MDF/g-AC
t2=A(:,1)/7.03;y2=A(:,2);%横坐标由g-MDF/g-AC化为h-1
q10=zeros(1,20);
c10=zeros(1,20);
q20=zeros(1,20);
c20=zeros(1,20);
y0=[c10 q10 c20 q20];
c0=7.03;
options=odeset('relTol',1e-6);
tspan=[0:0.02:t2(end)];
[t,y]=ode15s(@fangcheng,tspan,y0,options);
plot(t*7.03,y(:,40)/c0,'ro',t*7.03,y(:,80)/c0,'bo','LineWidth',3);hold on
plot(A(:,1),A(:,2),A(:,1),A(:,3),'LineWidth',3)
xlabel('Amount of treated MDF(g-MDF/g-AC)'),ylabel('c/c0'),axis([-0.05 120 -0.05 1.05]),
legend('拟合DBT','拟合MNA','实验曲线DBT','实验曲线MNA','Location','best') ,title('DBT+MNA-RS')
grid,hold off
%---------------------------------------------------------
function dydt=fangcheng(t,y)
kf=200;
a=0.363;%空隙率
u=0.361032/a;%床层颗粒间隙中流体流速
c0=7.03;
L=8.31e-2;
N=20;
dz=L/N;
p=436.8/(1-0.363);%颗粒表观密度对曲线影响较大从434.2该为436.8
q1=y(1:N);c1=y(N+1:2*N);
q2=y(2*N+1:3*N);c2=y(3*N+1:4*N);
[ce1,ce2]=equilibrium(q1,q2);
%-----------定义dq/dt-----------------------------
for j=1:N
dq1dt(j)=kf/p*(c1(j)-ce1(j));
dq2dt(j)=kf/p*(c2(j)-ce2(j));
end
%---------定义dC/dt-----------------------------------
dc1dt(1)=-u*(c1(1)-c0)/dz-(1-a)/a*p*dq1dt(1);
dc2dt(1)=-u*(c2(1)-c0)/dz-(1-a)/a*p*dq2dt(1);
for j=2:N
dc1dt(j)=-u*(c1(j)-c1(j-1))/dz-(1-a)/a*p*dq1dt(j);
dc2dt(j)=-u*(c2(j)-c2(j-1))/dz-(1-a)/a*p*dq2dt(j);
end
dydt=[dq1dt dc1dt dq2dt dc2dt]';
function [ce1,ce2]=equilibrium(q1,q2)
[ce1,ce2]=solve('q1=0.45967ce1^0.3839/(ce1^0.0457+0.3*ce2^0.0003206',...
'q2-0.17356*ce2^(-1.55201)/(ce2^(-2)+0.1*ce1^(-1))');[ Last edited by 小木虫zb on 2013-4-26 at 09:30 ]

回复此楼