数学实验6

实验6-非线性方程求解实验目的1.掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作初步分析。2.练习用非线性方程和方程组建立实际问题的模型并进行求解。实验内容Q6给定4种物质对应的参数ai,bi,ci和交互作用矩阵Q如下：a1=18.607a2=15.841a3=20.443a4=19.293b1=2643.31b2=2755.64b3=4628.96b4=4117.07c1=239.73c2=219.16c3=252.64c4=227.44Q=1.0000.1922.1691.6110.3161.0000.4770.5240.3770.3601.0000.2960.5240.2822.0651.000在压强p=760mmHg下，为了形成均相共沸混合物，温度和组分分别是多少？请尽量找出所有可能的解。【问题分析】设组分i所占比例为xi，则有:∑组分i饱和气相压强为pi，活度系数为有，据饱和气相压强与活度系数的计算式(设qij为组分i对组分j的交互作用系数)：∑∑(∑)将上述式子整合为：(∑)∑∑将式子变形为：((∑)∑(∑))【运算程序】首先编写如下的函数M文件：functionf=azeofun(XT,n,P,a,b,c,Q)x(n)=1;fori=1:n-1x(i)=XT(i);x(n)=x(n)-x(i);endT=XT(n);p=log(P);fori=1:nd(i)=x*Q(i,1:n)';dd(i)=x(i)/d(i);endfori=1:nf(i)=x(i)*(b(i)/(T+c(i))+log(x*Q(i,1:n)')+dd*Q(1:n,i)-a(i)-1+p);end然后用所给数据编程，作如下计算：n=4;P=760;a=[18.607,15.841,20.443,19.293]';b=[2643.31,2755.64,4628.96,4117.07]';c=[239.73,219.16,252.64,227.44]';Q=[1.00.1922.1691.6110.3161.00.4770.5240.3770.3601.00.2960.5240.2822.0651.0];XT0=[0.25,0.25,0.25,50];[XT,Y]=fsolve(@azeofun,XT0,[],n,P,a,b,c,Q)得到XT=[0.00000.58580.414271.9657]Y=1.0e-06*[-0.0009-0.04230.4428-0.4700]即4种物质组成均相共沸混合物时的比例分别为0,58.58%，41.42%，0，温度为71.9657°C。改变温度进行多次尝试，可以得到如下结果初值解XT0234T[0.25,0.25,0.25,50]0.00000.58580.41420.000071.9657[0.25,0.4,0.1,70]0.00000.78030.00000.219776.9613[0.5,0,0.5,100]0.00000.00000.00001.000097.7712【结果分析】本题与课本例题非常相似，原理及运算程序均比较简单。在做本题的过程中遇到最大的bug是设立不同的初值有可能得到相同的结果，所以在找不同的结果设立初值时花费了很长的时间。Q8假设商品在t时刻的市场价格为p（t），需求函数为D(p(t))=c-dp(t)（c，d≥0）。而生产方的期望价格为q（t），供应函数为S（q（t））。当供销平衡时S(q(t))=D（q(t)）。若期望价格与市场价格不符，商品市场不均衡，生产方t+1时期的期望价格将会调整，方式为q(t+1)-q(t)=r[p(t)-q(t)](0r1)，以p(t)=[c-D(p(t))]/d=[c-S(q(t))]/d带入,得到关于q(t)的递推方程.设S(x)=arctan(μx),μ=4.8,d=0.25,r=0.3,以c为可变参数,讨论期望价格q(t)的变化规律,是否有混沌现象出现,并找出前几个分岔点,观察分岔点的极限趋势是否符合Feigenbaum常数揭示的规律。【问题分析】依题意可得：q(t+1)=(1-r)q(t)+r(c-arctan(µq(t)))/d;代入数据得：q(t+1)=0.7q(t)+1.2(c-arctan(4.8q(t)));由此可构建函数并寻找规律。【运算程序】首先编写如下运算程序：functionchaos(iter_fun,x0,r,n)kr=0;forrr=r(1):r(3):r(2)kr=kr+1;y(kr,1)=feval(iter_fun,x0,rr);fori=2:n(2)y(kr,i)=feval(iter_fun,y(kr,i-1),rr);endendplot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');在本题中，迭代函数为：functiony=iter01(x,c)u=4.8;d=0.25;r=0.3;y=(1-r)*x-r*atan(u*x)/d+r*c/d;输入如下命令chaos(@iter01,5,[0.2,1.2,0.001],[100,200])可以得到如下图所示的分叉和混沌图（横坐标表示变化的参数，纵坐标表示迭代序列的极限）。从图中可以看出，分岔点的坐标为：1.0760.9400.9060.8952323243不难发现，n值越大，比值越接近Feigenbaum常数。再编写如下程序：functiony=iter02(c)u=4.8;d=0.25;r=0.3;q(1)=0.5;forn=1:1:50q(n+1)=(1-r)*q(n)-r*atan(u*q(n))/d+r*c/d;endN=1:1:51;plot(N,q);若考虑c的影响，只需将c的不同取值插入到程序中运行即可。1.c=1.2,有一个收敛子列，得到的图像如下：2.c=0.98时，有两个收敛子列。得到的图像如下：3．c=0.93时，有4个收敛子列。得到的图像如下：4.c=0.9时，有8个收敛子列。得到的图像如下：5.c=0.7时，处于混沌状态。得到的图像如下：【结果分析】本题可以用解析解法来解，进而验证。1.当方程有一个收敛极限时：根据q(t)=0.7q(t)+1.2(c-arctan(4.8q(t)))解得：q1.0769当方程有两个收敛极限时：根据q(t+2)=q(t)解得：1.0769q0.9065当方程有四个收敛极限时：根据q(t+4)=q(t)解得：0.9065q0.8968当方程有八个收敛极限时：根据q(t+8)=q(t)解得：0.8968q0.8685对比计算得到的结果和读图得到的结果会发现：1.读图难免有误差存在2.对特殊的c值单独作图的情况也与解析解的结果符合得很好。3.当n越大时，比值越接近于Feigenbaum常数。实验总结在做第一题的时候，一定要注意的地方就是迭代初值的选取。在具体的实践过程中，通过简单的画图方法来确定初始值的选取是一个解决问题的方法。经过适当的尝试，就可以找到一个落在收敛域内的迭代初值。在做第二题的时候，取c值的时候我依据的是读图的结果，那个时候得到的图像就有一些不满足要求，所以应该从取c值的时候就意识到读图有误差。后来用解析法检验的时候发现果然如此。我认为这次的MATLAB练习非常具有实用性，因为以前就解过共沸物组成的有关问题，以及解过收敛子列的题目。所以真的不得不说，MATLAB好强大，很复杂的过程都可以被轻松解答。通过这次课堂上的学习，我了解到借助计算机强大的计算功能，我们可以通过一定的数值分析，运用诸如牛顿法、拟牛顿法等方法，得到理想的近似解。\(^o^)/~