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

T=XT(n); p=log(P); for i=1:n d(i) = x * Q(i,1:n)'; dd(i)=x(i)/d(i); end for 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;
xn 1 xi
i 1
n 1
将它带入（2）式，得到含有 n 个未知数 XT=（ x1 , x2 ,..., xn1 , T ）的非线性方程组。 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(n)-x(i); end
，其中
是组分 i
的饱和汽相压强， 与温度 T 有关， 由 常数。 是组分 i 的液相活度系数，由
确定， 其中
为
确定，其中 qij 表示组分 i 与组分 j 的交互作用参数，qij 构成交互作用矩阵 Q，Q 不一定是对称矩阵。从而经过对数变换，稳定 条件可以改写成
（2） 则给定组成均相共沸混合物的 n 种物质，参数 ai, bi, ci 和交互作用矩阵 Q 是可以通过 实验得到的，可以作为已知系数。在一定的压强 P 下，模型（1），(2) 描述了确定均相共 沸混合物的组分 xi 的条件。 且（1），（2）是含有 n 1 个未知数非线性方程组，可以利用 MATLAB 优化工具包的 fsolve 求解，注意到（1）式是一个简单的线性等式，可以从中消去 1 个未知数，这通常会 使解的效果更好。例如，我们可以从（1）中解出
y(kc,i)=feval(ex_fun,y(kc,i-1),cc); end
end plot([c(1):c(3):c(2)],y(:,n(1)+1:n(2)),'k.') hold on fplot('0',[0,2],'r'); 在本题中，迭代函数是： function y=ex01(x,c) y=0.3*((c-atan(4.8*x))/0.25-x)+x; end 于是输入命令：chaosex(@ex01,0,[0,1.5,0.001],[100,200]) 得到 %做出 0 水平线
1 若������������ ，������������ 已知，且给定收获量ℎ������ 。根据 x
果：
b x
k 1
n
k
k 1 sk xk hk 可以得到以下结 和x
பைடு நூலகம்
1 x1 b1x1 b2 x2 ... bn1xn1 bn xn x 2 x2 s1 x1 h1 x 3 x3 s2 x2 h2 x
x1 -0.0000 -0.0000 -0.0000 -0.0000
经过尝试在 xi 0(i 1, 2,3, 4) 的条件下
当 XT0=[1/3,1/3,0,76]即假设没有第三种物质时，且初始温度设为 76℃时得到合理 解如下： XT =[-0.0000 0.7803 0.0000 0.0948 76.9613] -0.7068 0.2076]
3、问题总结：这道题目中解不是唯一的，这跟初始值的给定有关，而且还要注意每个组分 的含量是大于 0 的。所以对于初值要进行多次的尝试才能得到最终的答案。
二、混沌现象 1、问题分析：根据题目可以总结出期望价格的迭代式：
q(t 1) (1 r )q(t ) r (c arctan( q(t ))) / d
数学实验
郭明钊
实验 5 一、种群繁殖 问题分析：这道题和课本上的例题比较像，不同点在于多了收获量ℎ������ ，还有限制条件，即各
2012011880
化 21
k xk (K=1,2…,n)。 个年龄的种群数量保持不变， x 这道题的思路就是像例题那样列出线性代
数方程组然后解之。 1、 模型的建立：
获量都为 500.但是从结果看出，x5 为-259.49，但种群数量不可能为负数，在本题所给条件 下，无法使 h1~h5 均为 500。 4、 问题总结： 这道题用 matlab 实现时十分的容易， 就是解基本的线性代数方程组， 这道题的关键就是 如何正确地建立模型，只要找准数量之间的关系就能很好的建立模型。
观察上图，可用 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)'; end subplot(3,2,1),plot(k,xx(:,1)),gtext('c=0.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')得到图形如下：
代入数值后得：
q(t 1) 0.7q(t ) 1.2(c arctan( q(t )))
题目要求是观察是否有混沌现象并找出前几个分岔点，而且要观察是否符合 Feigenbaum 常数体现的规律。于是就应该先做出混沌分岔图 2、Matlab 做出混沌分岔图，首先编写如下程序： function chaosex(ex_fun,x0,c,n) x0 是迭代初值； kc=0; for cc=c(1):c(3):c(2) 长 kc=kc+1; y(kc,1)=feval(ex_fun,x0,cc); for i=2:n(2) 值被舍弃 %输入中 n(2)是迭代序列的长度，但画图时前 n(1)个迭代 % 输入中[r(1),r(2)]是参数变化的范围，r(3) 是步 % 该函数没有返回值； iter_fun 是迭代函数 （句柄） ；
Y =1.0e-010 *[ -0.0004
也就是 4 种物质含量分别是 0.00%，58.58%，41.42%，0.00%温度是 71.9657℃ 在上述计算中初值 XY 0 的取法是：四种物质各占约 1/4,温度为 70℃。如果取
其他值会如何呢？
初值 XT0 [0,1/3,1/3,76] [1/3,0,1/3,76] [1/3,1/3,0,76] [1/3,1/3,1/3,76] 解 x2 x3 x4 T 0.5858 0.4124 0.0000 71.9657 -0.0000 1.0000 0.0000 82.5567 0.7803 0.0000 0.2197 76.9613 0.5858 0.4142 0.0000 71.9656
Y =1.0e-006 *[0.0000
此时混合物中包含有 78.03%的第二种物质和 21.97%的第四种物质。 温度为 76.9613℃ 综上所述，此题目中所有的解为
1 4 种物质含量分别是 0.00%，58.58%，41.42%，0.00%温度是 71.9657℃ ○
2 4 种物质含量分别是 0.00%，78.03%，0.00%，21.97%温度是 76.9613℃ ○
bn 1 0 0 1
sn 1
bn 0 0 0 1
当年各个年龄段的数量（向量）为： X [ x1 , x2 ,..., xn1, xn ]T 再令 H [0, h1 , h2 ,..., hn2 , hn1 ]T 于是就有： H LX 2、 对于第二问，运用 matlab 实现： clc L=[-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; 400; 200; 100]; x=L\H; x1=vpa(x,5) 得到结果如下： %求解当前时段各年龄组的数量，保留 5 位有效数字。 x1 = 8481.0 2892.4 1335.4 601.27 140.51 于是可以知道 x1 ~ x5 分别是 8481.0 2892.4 1335.4 601.3 140.5 3、 对于第三问，只和第二问有一点差别，就是 h1 ~ h5 均是 500.用 matlab 实现时和第二问 的方法一样。 clc L=[-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=L\H; x1=vpa(x,5) %求解当前时段各年龄组的数量，保留 5 位有效数字。 得到的结果如下： x1 = 10981.0 3892.4 1835.4 601.27 -259.49 于是 x1 ~ x5 分别是 10981.0 3892.4 1835.4 601.3 -259.5 时可以实现各个年龄段的收
实验 6 一、均相共沸混合物的组分 1、问题分析：所谓共沸混合物，是指有两种或两种以上物质组成的液体混合物，在某 种压力下被蒸馏或局部汽化时， 在气体状态下和在液体状态下保持相同的组分。 课本中 已经建立了模型如下： 组分 xi（i=1，…, n）应满足：
（1）
在压强 P 不大的情况下，稳定条件可以表示
……
n xn sn1xn1 hn1 x
稍作变换可得：
(b1 1) x1 b2 x2 ... bn1xn1 bn xn 0 s1 x1 x2 h1 s2 x2 x3 h2
……
sn1 xn1 xn hn1
于是可以得到模型如下：
b1 1 b2 b3 1 0 s1 0 s2 1 L= 0 0 0 0 0 0
由此可见 c=1.2 时有一个收敛极限，c=1.0 时有两个收敛极限，c=0.92 时有四个收敛极 限，当 c=0.9 时，有八个收敛极限，c=0.88 和 0.8 时看起来就比较混乱了，可以认为出现了 混沌现象。 下面分析以上 c 为不同值时的分岔现象及具体分岔点： 这些点在分岔混沌图上面通过定点显示坐标的方式就能够找出来。 前几个分岔点是：1.0769,0.949, 0.9065,0.8968，0.8685 当然这也能够通过代数运算的方式计算出来，理论上有以下结论： （1） 方程只有一个收敛极限时： 根据 q 0.7q 1.2(c arctan( q)) 求得 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 但是由于计算分析的能力有限以上的理论计算没能实现。
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.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 -0.0446 71.9657] 0.2627 -0.1489]