luruiliu
2007-12-26, 11:10
非常感谢你能打开这个帖子，我是matlab的初学者，以前在学校学过用过一下子matlab，现在基本上还给我的老师了，工作遇到了要用matlab，有些问题不懂，想让高手指点一下子，就是matlab中的小程序如何衔接调用，譬如附件中的几个matlab的程序，基本上能看懂，但是就是不知道怎么衔接起来，还有一些设定参数不知道怎么改变调整。清高手指点一下子，见笑了。谢谢。
1.稳态解程序
1.1函数wt.m〔稳态典型参数的输入)
global YHY AF lIX B
global UH KSB H
global UAK NH KOAB A
global KH
global S10S SOX IOX SOX BHOX BAO0 0S NOOS NHOS ALKO %进水组分
global VQ dh dx SRT O Os
%模型参数值
YH =0.71;
YA =0.24;
FI= 0.07;
IX B=0.086;
UH=5.70;
KS=20.0;
BH=0.48;
UA=0.45;
KNH=0.5;
KOA=0.3;
BA=o.1;
KH=100;
%进水划分
S10 =40;
SS O=64;
X1 0=40;
XS O=160;
XB HO=96;
XB AO=0.001;
SN OO=1;
00 =0;
SN HO=12.5;
SA LKO=6.0;
%运行参数
V= 45 ;
Q= 10 0;
SR T= 10 ;
Os =8 .6 37;
0= 2;
dx =工/SR丁;
dh =V /Q ;
1.2在运行wt二的基础上运行chengxu.m 数，chengxu.m如下:
global YH YA FI IXB
global UH KS BH
global UA KNH KOA BA
globalK H
global S11 SSl XII XSl XBHI XBA1 Ol SNOI SNHI SALK1
global V Q dh dx Os O
syms SI SS XI XS XBH XBA SNO SNH SALK Kla
component0=[SISS X IX SX BHX BAK laS NOS NHS ALK];
influent=[SI] SS1 XII XSl XBHl XBAI Ol SNOI SNHI SALKI];
velocity=sym(zeros(1,10));
table=sym(zeros(5,10));
process=sym(zeros(5,1));
equation=sym(zeros(10,1));
process(1)=UH*(S S/(KS+SS))*XBH;
process(2)= UA*(SNH/(KNH+SNH))*(O/(KOA+O))*XBA;
process(3)=BH*XBH;
process(4)=BA*XBA;
process(5)=KH* XS;
table(1,2)=-1/Y H ;
table(1,5)=I;
table( 1,7)=-(1-YH)/YH;
table(1,9)=-IXB;
table(1,10)=-IXB/14;
table(2,6)=1;
table(2,7)=-(4.57-YA)/YA;
table(2,8)=1/YA;
talbe(2,9)=-IXB-1/YA;
table(2,10)=-IXB /14 -1/(7*YA) ;
table(3,3)=FI;
table(3,4)=1-FI;
table(3,5)=-1;
talbe(3,9)=IXB;
table(4,3)=FI;
table(4,4)=1-FI;
table(4,6)=l;
table(4,9)=IXB;
table(5,2)=1;
table(5,4)=-1;
for i=1:10
for j= 1:5
velocity(i)=velocity(i)+table(j,i)*process(j);
end
end
v=velocity';
equation(1,1)=influent(1)*dh-component0(1)*dh+velocity(1);
equation(2,1)=influent(2)*dh-component0(2)*d h+velocity(2);
equation(3,1)=influent(3)*dh-component0(3)*dx+velocity(3);
equation(4,1)=influent(4)*dh-component0(4)*dx+velocity(4);
equation(5,1)=influent(5)*dh-component0(5)*dx+velocity(5);
equation(6,1)=influent(6)*dh-component0(6)*dx+velocity(6);
equation(7,1)=influent(7)*dh+component0(7)*(Os-O)+velocity(7);
equation(8,1)=influent(8)*dh-component0(8)*dh+velocity(8);
equation(9,1)=influent(9)*dh-component0(9)*dh+velocity(9);
equation(10,1)=influent(l0)*dh-component0(10)*dh+velocity(10);
equation;
Sl='x(1)';
SS='x(2)';
XI='x(3)';
XS='x(4)';
XBH='x(5)';
XBA='x(6)';
Kla='x(7)';
SNO='x(8)';
SNH='x(9)';
SALK='x(10)';
subs(equation(l,l))
subs(equation(2,1))
subs(equation(3,1))
subs(equation(4,1))
subs(equation(5,1))
subs(equation(6,1))
subs(equation(7,1))
subs(equation(8,1))
subs(equation(9,1))
subs(equation(10,1))
1.3运行chengxu.m后，得到10个方程，组成非线性常微分方程组，用i'数rigid.m表示:
%rigid.m
function dx=rigid(t,x)
dx=zeros(10,1)
dx(1)=800/9-20/9*x(1);
dx(2)=1280/9-20/9*x(2)-570/71*x (2)/(20+x(2))*x(5)+100*x(4);
dx(3)=800/9-1/10*x(3)+21/625*x(5)+7/1000*x(6);
dx(4)=3200/9-1001/10*x(4)+279/625*x(5)+93/1000*x(6);
dx(5)=640/3-29/50*x(5)+57/10*x(2)/(20+x(2))*x(5);
dx(6)=1/4 50-1/5 *x (6)+9/23* x( 9)/(1/2+x(9))*x( 6);
dx(7)=6637/10 00*x(7)-1653/710*x(2)/(20+x(2))*x(5)-1299/184*x(9)/(1/2 +x(9))*x(6);
dx(8)=-20/9*x(8)+75/46*x(9)/(1/2 +x(9))*x(6);
dx(9)=250/9-20/9*x(9)-2451/5 000*x(2)/(20+x(2))*x(5)+43/5000*x(6);
dx(10)=40/3-20/9*x(10)-2451170000*x(2)/(20+x(2))*x(5)-37887/161000*x(9)/(1/2+x(9))*x(6);
1.4 用龙格一库塔法解上面的常微分方程组，
%main.m
xO=[40 1.5 889 30 1450 90 150 40
[T,X]=ode45('rigid',[O 0.001],x0)
程序为main.m:
1.5 程序运行结果如下
T=
0
0.0 0 00
0.0 0 0 1
(省略)
0.0009
0.0009
0.0010
0.0010
0.0010
Oe+003*
0.0 4 00 0.0015
0.0003 0.0025
0.8890 0.0300 1.4500 0 .0900 0 .1500 0 .0400
0.0 4 00
0.0003
0.0 0 16
0.0025
0.8890 0.0300 4500 0 .0900 0 .1500 0 .0400
0.0 4 00
0.0003
0.0 0 16
0.0025
0.8890 0.0299 4500 0 .0900 0 .1500 0 .0400
0.0 4 00 0.0032
0.0003 0.0025
0.0 4 00 0 .0032
0.0003 0.0025
0.0 40 0 0 .0032
0.0003 0.0025
0.0 4 00 0 .0033
0.0003 0.0025
(省略)
0.8890 0.0282 1.4502 0 .0900 0 .1504 0 .0400
0.8890 0.0282 1.4503 0.0900 0.1504 0 .0400
0.8890 0.0281 1.4503 0.0900 0.1504 0 .0400
0.8890 0.0281 1.4503 0 .0900 0 .1504 0 .0400
118
1.稳态解程序
1.1函数wt.m〔稳态典型参数的输入)
global YHY AF lIX B
global UH KSB H
global UAK NH KOAB A
global KH
global S10S SOX IOX SOX BHOX BAO0 0S NOOS NHOS ALKO %进水组分
global VQ dh dx SRT O Os
%模型参数值
YH =0.71;
YA =0.24;
FI= 0.07;
IX B=0.086;
UH=5.70;
KS=20.0;
BH=0.48;
UA=0.45;
KNH=0.5;
KOA=0.3;
BA=o.1;
KH=100;
%进水划分
S10 =40;
SS O=64;
X1 0=40;
XS O=160;
XB HO=96;
XB AO=0.001;
SN OO=1;
00 =0;
SN HO=12.5;
SA LKO=6.0;
%运行参数
V= 45 ;
Q= 10 0;
SR T= 10 ;
Os =8 .6 37;
0= 2;
dx =工/SR丁;
dh =V /Q ;
1.2在运行wt二的基础上运行chengxu.m 数，chengxu.m如下:
global YH YA FI IXB
global UH KS BH
global UA KNH KOA BA
globalK H
global S11 SSl XII XSl XBHI XBA1 Ol SNOI SNHI SALK1
global V Q dh dx Os O
syms SI SS XI XS XBH XBA SNO SNH SALK Kla
component0=[SISS X IX SX BHX BAK laS NOS NHS ALK];
influent=[SI] SS1 XII XSl XBHl XBAI Ol SNOI SNHI SALKI];
velocity=sym(zeros(1,10));
table=sym(zeros(5,10));
process=sym(zeros(5,1));
equation=sym(zeros(10,1));
process(1)=UH*(S S/(KS+SS))*XBH;
process(2)= UA*(SNH/(KNH+SNH))*(O/(KOA+O))*XBA;
process(3)=BH*XBH;
process(4)=BA*XBA;
process(5)=KH* XS;
table(1,2)=-1/Y H ;
table(1,5)=I;
table( 1,7)=-(1-YH)/YH;
table(1,9)=-IXB;
table(1,10)=-IXB/14;
table(2,6)=1;
table(2,7)=-(4.57-YA)/YA;
table(2,8)=1/YA;
talbe(2,9)=-IXB-1/YA;
table(2,10)=-IXB /14 -1/(7*YA) ;
table(3,3)=FI;
table(3,4)=1-FI;
table(3,5)=-1;
talbe(3,9)=IXB;
table(4,3)=FI;
table(4,4)=1-FI;
table(4,6)=l;
table(4,9)=IXB;
table(5,2)=1;
table(5,4)=-1;
for i=1:10
for j= 1:5
velocity(i)=velocity(i)+table(j,i)*process(j);
end
end
v=velocity';
equation(1,1)=influent(1)*dh-component0(1)*dh+velocity(1);
equation(2,1)=influent(2)*dh-component0(2)*d h+velocity(2);
equation(3,1)=influent(3)*dh-component0(3)*dx+velocity(3);
equation(4,1)=influent(4)*dh-component0(4)*dx+velocity(4);
equation(5,1)=influent(5)*dh-component0(5)*dx+velocity(5);
equation(6,1)=influent(6)*dh-component0(6)*dx+velocity(6);
equation(7,1)=influent(7)*dh+component0(7)*(Os-O)+velocity(7);
equation(8,1)=influent(8)*dh-component0(8)*dh+velocity(8);
equation(9,1)=influent(9)*dh-component0(9)*dh+velocity(9);
equation(10,1)=influent(l0)*dh-component0(10)*dh+velocity(10);
equation;
Sl='x(1)';
SS='x(2)';
XI='x(3)';
XS='x(4)';
XBH='x(5)';
XBA='x(6)';
Kla='x(7)';
SNO='x(8)';
SNH='x(9)';
SALK='x(10)';
subs(equation(l,l))
subs(equation(2,1))
subs(equation(3,1))
subs(equation(4,1))
subs(equation(5,1))
subs(equation(6,1))
subs(equation(7,1))
subs(equation(8,1))
subs(equation(9,1))
subs(equation(10,1))
1.3运行chengxu.m后，得到10个方程，组成非线性常微分方程组，用i'数rigid.m表示:
%rigid.m
function dx=rigid(t,x)
dx=zeros(10,1)
dx(1)=800/9-20/9*x(1);
dx(2)=1280/9-20/9*x(2)-570/71*x (2)/(20+x(2))*x(5)+100*x(4);
dx(3)=800/9-1/10*x(3)+21/625*x(5)+7/1000*x(6);
dx(4)=3200/9-1001/10*x(4)+279/625*x(5)+93/1000*x(6);
dx(5)=640/3-29/50*x(5)+57/10*x(2)/(20+x(2))*x(5);
dx(6)=1/4 50-1/5 *x (6)+9/23* x( 9)/(1/2+x(9))*x( 6);
dx(7)=6637/10 00*x(7)-1653/710*x(2)/(20+x(2))*x(5)-1299/184*x(9)/(1/2 +x(9))*x(6);
dx(8)=-20/9*x(8)+75/46*x(9)/(1/2 +x(9))*x(6);
dx(9)=250/9-20/9*x(9)-2451/5 000*x(2)/(20+x(2))*x(5)+43/5000*x(6);
dx(10)=40/3-20/9*x(10)-2451170000*x(2)/(20+x(2))*x(5)-37887/161000*x(9)/(1/2+x(9))*x(6);
1.4 用龙格一库塔法解上面的常微分方程组，
%main.m
xO=[40 1.5 889 30 1450 90 150 40
[T,X]=ode45('rigid',[O 0.001],x0)
程序为main.m:
1.5 程序运行结果如下
T=
0
0.0 0 00
0.0 0 0 1
(省略)
0.0009
0.0009
0.0010
0.0010
0.0010
Oe+003*
0.0 4 00 0.0015
0.0003 0.0025
0.8890 0.0300 1.4500 0 .0900 0 .1500 0 .0400
0.0 4 00
0.0003
0.0 0 16
0.0025
0.8890 0.0300 4500 0 .0900 0 .1500 0 .0400
0.0 4 00
0.0003
0.0 0 16
0.0025
0.8890 0.0299 4500 0 .0900 0 .1500 0 .0400
0.0 4 00 0.0032
0.0003 0.0025
0.0 4 00 0 .0032
0.0003 0.0025
0.0 40 0 0 .0032
0.0003 0.0025
0.0 4 00 0 .0033
0.0003 0.0025
(省略)
0.8890 0.0282 1.4502 0 .0900 0 .1504 0 .0400
0.8890 0.0282 1.4503 0.0900 0.1504 0 .0400
0.8890 0.0281 1.4503 0.0900 0.1504 0 .0400
0.8890 0.0281 1.4503 0 .0900 0 .1504 0 .0400
118