实验5,6 线性代数方程组的数值解法非线性方程组求解

1、数学实验郭明钊 化21实验5一、 种群繁殖问题分析：这道题和课本上的例题比较像，不同点在于多了收获量hk，还有限制条件，即各个年龄的种群数量保持不变，(K=1,2,n)。这道题的思路就是像例题那样列出线性代数方程组然后解之。1、 模型的建立：若bk，sk已知，且给定收获量hk。根据和可以得到以下结果： 稍作变换可得：于是可以得到模型如下： L=当年各个年龄段的数量（向量）为：再令 于是就有： 2、 对于第二问，运用matlab实现：clcL=-1 0 5 3 0;0.4 -1 0 0 0;0 0.6 -1 0 0;0 0 0.6 -1 0;0 0 0 0.4 -1;%构造L矩阵H=0; 500

2、; 400; 200; 100;x=LH;x1=vpa(x,5) %求解当前时段各年龄组的数量，保留5位有效数字。 得到结果如下：x1 = 8481.0 2892.4 1335.4 601.27 140.51 于是可以知道分别是8481.0 2892.4 1335.4 601.3 140.53、 对于第三问，只和第二问有一点差别，就是均是500.用matlab实现时和第二问的方法一样。clcL=-1 0 5 3 0;0.4 -1 0 0 0;0 0.6 -1 0 0;0 0 0.6 -1 0;0 0 0 0.4 -1;%构造L矩阵H=0; 500; 500; 500; 500;x=LH;x1=

3、vpa(x,5) %求解当前时段各年龄组的数量，保留5位有效数字。得到的结果如下：x1 = 10981.0 3892.4 1835.4 601.27 -259.49于是分别是10981.0 3892.4 1835.4 601.3 -259.5 时可以实现各个年龄段的收获量都为500. 但是从结果看出，x5为-259.49，但种群数量不可能为负数，在本题所给条件下，无法使h1h5均为500。4、 问题总结：这道题用matlab实现时十分的容易，就是解基本的线性代数方程组，这道题的关键就是如何正确地建立模型，只要找准数量之间的关系就能很好的建立模型。实验6一、 均相共沸混合物的组分1、问题分析：所

4、谓共沸混合物，是指有两种或两种以上物质组成的液体混合物，在某种压力下被蒸馏或局部汽化时，在气体状态下和在液体状态下保持相同的组分。课本中已经建立了模型如下：组分xi（i=1，, n）应满足： （1）在压强P不大的情况下，稳定条件可以表示，其中是组分i的饱和汽相压强，与温度T有关，由确定，其中为常数。是组分i的液相活度系数，由 确定，其中qij表示组分i与组分j的交互作用参数，qij构成交互作用矩阵Q，Q不一定是对称矩阵。从而经过对数变换，稳定条件可以改写成 （2）则给定组成均相共沸混合物的n种物质，参数ai, bi, ci和交互作用矩阵Q是可以通过实验得到的，可以

5、作为已知系数。在一定的压强P下，模型（1），(2) 描述了确定均相共沸混合物的组分xi的条件。且（1），（2）是含有个未知数非线性方程组，可以利用MATLAB优化工具包的fsolve求解，注意到（1）式是一个简单的线性等式，可以从中消去1个未知数，这通常会使解的效果更好。例如，我们可以从（1）中解出将它带入（2）式，得到含有n个未知数XT=（）的非线性方程组。2、matlab实现，根据题目中所给的数据可以matlab编写如下：首先编写如下的函数M文件：function f=azeofun(XT,n,P,a,b,c,Q)x(n)=1;for i=1:n-1 x(i)=XT(i); x(n)=x(

6、n)-x(i);endT=XT(n);p=log(P);for i=1:n d(i) = x * Q(i,1:n)' dd(i)=x(i)/d(i);endfor i=1:n f(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&

7、#39;Q=1.0 0.192 2.169 1.611 0.316 1.0 0.477 0.524 0.377 0.360 1.0 0.296 0.524 0.282 2.065 1.0;XT0=1/4,1/4,1/4,70;XT,Y=fsolve(azeofun,XT0,n,P,a,b,c,Q)得到XT =0.0000 0.5858 0.4142 71.9657Y =1.0e-010 * -0.0004 -0.0446 0.2627 -0.1489也就是4种物质含量分别是0.00%，58.58%，41.42%，0.00%温度是71.9657在上述计算中初值的取法是：四种物质各占约1/4,温度

8、为70。如果取其他值会如何呢？经过尝试在的条件下当XT0=1/3,1/3,0,76即假设没有第三种物质时，且初始温度设为76时得到合理解如下：XT =-0.0000 0.7803 0.0000 76.9613Y =1.0e-006 *0.0000 0.0948 -0.7068 0.2076此时混合物中包含有78.03%的第二种物质和21.97%的第四种物质。温度为76.9613综上所述，此题目中所有的解为 4种物质含量分别是0.00%，58.58%，41.42%，0.00%温度是71.9657 4种物质含量分别是0.00%，78.03%，0.00%，21.97%温度是76.96133、问题总结

9、：这道题目中解不是唯一的，这跟初始值的给定有关，而且还要注意每个组分的含量是大于0的。所以对于初值要进行多次的尝试才能得到最终的答案。二、混沌现象1、问题分析：根据题目可以总结出期望价格的迭代式：代入数值后得：题目要求是观察是否有混沌现象并找出前几个分岔点，而且要观察是否符合Feigenbaum 常数体现的规律。于是就应该先做出混沌分岔图2、Matlab做出混沌分岔图，首先编写如下程序：function chaosex(ex_fun,x0,c,n) % 该函数没有返回值；iter_fun是迭代函数（句柄）；x0是迭代初值；kc=0; for cc=c(1):c(3):c(2) % 输入中r(1

10、),r(2)是参数变化的范围，r(3) 是步长 kc=kc+1; y(kc,1)=feval(ex_fun,x0,cc); for i=2:n(2) %输入中n(2)是迭代序列的长度，但画图时前n(1)个迭代值被舍弃 y(kc,i)=feval(ex_fun,y(kc,i-1),cc); endendplot(c(1):c(3):c(2),y(:,n(1)+1:n(2),'k.') hold onfplot('0',0,2,'r'); %做出0水平线在本题中，迭代函数是：function y=ex01(x,c)y=0.3*(c-atan(4.8*

11、x)/0.25-x)+x;end于是输入命令：chaosex(ex01,0,0,1.5,0.001,100,200) 得到观察上图，可用c=1.2，1.0，0.92，0.90，0.88，0.80画图，matlab编写如下：c=0.8 0.88 0.9 0.92 1.0 1.2;x=0;n=40;for j=1:6 C=c(j); for i=1:n x(i+1)=0.3*(C-atan(4.8*x(i)/0.25-x(i)+x(i); end xx(:,j)=x' k=(0:40)'endsubplot(3,2,1),plot(k,xx(:,1),gtext('c=0.

12、8')subplot(3,2,2),plot(k,xx(:,2),gtext('c=0.88')subplot(3,2,3),plot(k,xx(:,3),gtext('c=0.9')subplot(3,2,4),plot(k,xx(:,4),gtext('c=0.92')subplot(3,2,5),plot(k,xx(:,5),gtext('c=1.0')subplot(3,2,6),plot(k,xx(:,6),gtext('c=1.2') 得到图形如下：由此可见c=1.2时有一个收敛极限，c=1.0

13、时有两个收敛极限，c=0.92时有四个收敛极限，当c=0.9时，有八个收敛极限，c=0.88和0.8时看起来就比较混乱了，可以认为出现了混沌现象。下面分析以上c为不同值时的分岔现象及具体分岔点：这些点在分岔混沌图上面通过定点显示坐标的方式就能够找出来。前几个分岔点是：1.0769, 0.949, 0.9065, 0.8968，0.8685当然这也能够通过代数运算的方式计算出来，理论上有以下结论：（1） 方程只有一个收敛极限时：根据求得c>1.0769（2） 方程只有两个收敛极限时：根据q（t+2）=q（t）可以解得：1.0769>c>0.949（3） 方程有四个收敛极限时：根据q（t+4）=q（t）可以解得：0.949>c>0.9065（4） 方程有八个收敛极限时：根据q（t+8）=q（t）可以解得：0.9065>c>0.8968但是由于计算分析的能力有限以上的理论