娃娃小囡 2008-8-21 14:33
帮帮小女啊，matlab解积分方程组
[b][font=宋体][size=5]四个方程联立求四个未知数，分别是：x1,x2,P1,P2,其余都是已知的
[/size][size=14pt]有三个方程非常简单，就是第一个左右两边都是积分，要求的未知数在积分上限里。牛人帮我解一下啊！！下面是我已经编出来的一部分，但是运行不对，出错，牛人只要帮我稍稍修改一下[/size][/font][/b]
[font=宋体][size=14pt]∫(0，P1)(n1/p)dp=[font=宋体][size=14pt]∫(0，P2)(n2/p)dp（其中n1，n2都可以用p表示：n1=n01*(b1*p^q1)/(1+*(b1*p)^q1),[/size][/font][/size][/font]
[font=宋体][size=14pt][font=宋体][size=14pt]n2=n02*(b2*p^q2)/(1+*(b2*p)^q2)）[/size][/font][/size][/font]
[font=宋体][size=14pt][font=宋体][size=14pt]P*y1=P1*x1[/size][/font][/size][/font]
[font=宋体][size=14pt][font=宋体][size=14pt]P*y2=P2*x2[/size][/font][/size][/font]
[font=宋体][size=14pt][font=宋体][size=14pt]x1+x2=1[/size][/font][/size][/font]
[font=宋体][size=14pt]
[/size][/font]
[font=宋体][size=14pt]%计算π*
%π*=∫(0，p)(nt/p)dp
%使用Matlab的符号积分
for i=1: L
aa=int('(b*p)^q/(1+(b*p)^q)*n0/p','p',0,'p0');
t=subs(aa,'b',m(2));
tt=subs(t,'q',m(3));
ttt=subs(tt,'n0',m(1));
pai(i)=subs(ttt,'p0',p(i));
pai(i)=double(pai(i));
end[/size][/font]
[font=宋体][size=14pt]%计算p1，p2，需要使用到纯组分的拟合表达式
n0=[5.0475 5.2755];
b=[18.4504 1.6382];
q=[0.6006 0.7519];%ethane, methane
for i=1: L-1
p1=1+0.3*i;
options=OPTIMSET('Tolx',1e-8,'Tolfun',1e8,'Display','iter','LargeScale','off','LevenbergMarquardt','on');
[p10(i),fval,flag]=fzero('fun2',p01,options,n0(2),b(2),q(2),pai(i+1));%methane
a3(i)=flag;
p2=0.3+0.05*i;
[p20(i),fval,flag]=fzero('fun2',p02,options,n0(1),b(1),q(1),pai(i+1));%ethane
b3(i)=flag;
end[/size][/font]
[font=宋体][size=14pt]function F=fun2(z,n0,b,q,I)[/size][/font]
[font=宋体][size=14pt]s1='(b*p)^q/(1+(b*p)^q)*n0/p';[/size][/font]
[font=宋体][size=14pt]aa(1)=int(s1,'p',0,'p10');[/size][/font]
[font=宋体][size=14pt]t(1)=subs(aa(1),'b',b);
tt(1)=subs(t(1),'q',q);
ttt(1)=subs(tt(1),'n0',n0);[/size][/font]
[font=宋体][size=14pt]a1=ttt(1);[/size][/font]
[font=宋体][size=14pt]a1=subs(a1,'p10',z);[/size][/font]
[font=宋体][size=14pt]a1=double(a1);[/size][/font]
[font=宋体][size=14pt]
F=I-a1;[/size][/font]