试验5线性代数方程组的数值解法
线性代数方程组的数值解法讨论
线性代数方程组的数值解法讨论解线性方程组的方法，主要分为直接方法和迭代方法两种。
直接法是在没有舍入误差的假设下能在预定的运算次数内求得精确解。
而实际上，原始数据的误差和运算的舍入误差是不可以避免的，实际上获得的也是近似解。
迭代法是构造一定的递推格式，产生逼近精确解的序列。
对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，因此比较受工程人员青睐。
小组成员本着工程应用，讨论将学习的理论知识转变为matlab 代码。
讨论的成果也以各种代码的形式在下面展现。
1 Jacobi 迭代法使用Jacobi 迭代法，首先必须给定初始值，其计算过程可以用以下步骤描述： 步骤1 输入系数矩阵A ，常熟向量b ，初值(0)x ，误差限ε，正整数N ，令1k =.步骤2 (0)11ni i ij jj ii j i x b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，(0)j x 代表(0)x 的第j 个分量。
步骤3 计算11ni i ij j j ii j i y b a x a =≠⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦∑，判断1max i i i n x y ε≤≤-<，如果是，则结束迭代，转入步骤5；否则，转入步骤4。
步骤4 判断k N =？如果是，则输出失败标志；否则，置1k k =+，i i x y ⇐，1,2,,i n =，转入步骤2。
步骤5 输出12,,n y y y 。
雅可比迭代代码function [x,k]=Fjacobi(A,b,x0,tol)% jacobi 迭代法 计算线性方程组% tol 为输入误差容限，x0为迭代初值max1= 300; %默认最多迭代300，超过要300次给出警告 D=diag(diag(A)); L=-tril(A,-1);U=-triu(A,1); B=D\(L+U); f=D\b; x=B*x0+f;k=1; %迭代次数while norm(x-x0)>=tol x0=x;x=B*x0+f; k=k+1;if(k>=max1)disp('迭代超过300次，方程组可能不收敛'); return; end%[k x'] %显示每一步迭代的结果 End2 高斯赛德尔迭代由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值，若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量(1)k i x +时，用最新分量11()k x +，12()k x +…(1)1k i x +-代替旧分量)1(k x ', )2(k x …)3(k x 就得到高斯赛德尔迭代格式，其数学表达式为：1(1)(1)()111(1,2,,)i n k k k ii ij j ij j j j i ii xb a x a x i n a -++==+⎛⎫=--= ⎪⎝⎭∑∑具体形式如下：()()()(1)()()()11221331111(1)(1)()()22112332222(1)(1)(1)(1)(1)112233,11111k k k k n n k k k k n n k k k k k n n n n n n n n nnx a x a x a x b a x a x a x a x b a x a x a x a x a x b a ++++++++--=----+=----+⋯⋯⋯⋯⋯⋯=-----+矩阵形式表示为：()(1)1(1)()(0,1,2,,),k k k k n +-+=++=x D Lx Ux b将(1)(1)()(0,1,2,,)k k k k n ++=++=Dx Lx Ux b 移项整理得： (1)1()1()()(0,1,2,,))k k x D L Ux D L b k n +--=-+-=记11(),()--=-=-M D L U g D L b ，则(1)()k k x x +=+M g高斯塞德尔迭代function [x,k]=Fgseid(A,b,x0,tol)%高斯-塞德尔迭代法 计算线性方程组 % tol 为误差容限max1= 300; %默认最高迭代300次D=diag(diag(A)); L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x=G*x0+f;k=1; while norm(x-x0)>=tol x0=x;x=G*x0+f; k=k+1;if(k>=max1)disp('迭代次数太多，可能不收敛'); return; end% [k,x'] %显示每一步迭代结果 End3 超松弛迭代法在工程中最常遇到的问题便是线性代数方程组的求解,而线性代数方程组的求解一般可以分为两类,一类是直接法(精确法),包括克莱姆法则方法、LD 分解法等，另一类是迭代法（近似法），包括雅克比迭代法、高斯迭代法、超松弛迭代法等。
实验5_线性代数方程组的数值解法
实验5 线性代数方程组的数值解法化工系 毕啸天 2010011811【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】题目3已知方程组Ax=b ，其中A ，定义为⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------=32/14/12/132/14/14/12/132/14/14/12/132/14/12/13A试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求：（1）选取不同的初始向量x (0)和不同的方程组右端项向量b ，给定迭代误差要求，用雅可比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论；（2）取定右端向量b 和初始向量x (0)，将A 的主对角线元素成倍增长若干次，非主对角线元素不变，每次用雅可比迭代法计算，要求迭代误差满足 ，比较收敛速度，分析现象并得出你的结论。
3.1 模型分析选取初始向量x(0) =(1,1,…,1)T ，b=(1,1,…,1)T ，迭代要求为误差满足 ，编写雅各比、高斯-赛德尔迭代法的函数，迭代求解。
3.2 程序代码function x = Jacobi( x0,A,b,m ) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B1=D\(L+U); f1=D\b; x(:,1)=x0;x(:,2)=B1*x(:,1)+f1; k=1;while norm((x(:,k+1)-x(:,k)),inf)>m x(:,k+2)=B1*x(:,k+1)+f1; k=k+1; endendfunction x = Gauss( x0,A,b,m )D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B2=(D-L)\U;f2=(D-L)\b;x(:,1)=x0;x(:,2)=B2*x(:,1)+f2;k=1;while norm((x(:,k+1)-x(:,k)),inf)>mx(:,k+2)=B2*x(:,k+1)+f2;k=k+1;endendA1=3.*eye(20,20);A2=sparse(1:19,2:20,-1/2,20,20);A3=sparse(1:18,3:20,-1/4,20,20);AA=A1+A2+A3+A2'+A3';A=full(AA);b=ones(20,1); %输入自选右端项向量b x0=ones(20,1); %输入自选初始向量x0 m=1e-5;x1=Jacob(x0,A,b,m);x2=Gauss(x0,A,b,m);结果输出数据：3.3.1将x0各分量初值置为0【分析】从数据中可以看出，当迭代的初值变化了，达到相同精度所需要的迭代次数也变化了。
线性方程组的数值解法实验报告
实验报告——线性方程组的数值解法姓名：班级：学号：日期：一 实践目的1. 熟悉求解线性方程组的有关理论和方法。
2. 会编列主元消去法，全主元消去法，雅克比迭代法及高斯-赛德尔迭代法的程序。
3.通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。
4. 进一步应用数学知识，扩展数学思维，提高编程能力。
二 问题定义及题目分析1.求解线性方程是实际中常遇到的问题。
而一般求解方法有两种，一个是直接求解法，一个是迭代法。
2.直接法是就是通过有限步四则运算求的方程准确解的方法。
但实际计算中必然存在舍入误差，因此这种方法只能得到近似解。
3.迭代法是先给一个解的初始近似值，然后按一定的法则求出更准确的解，即是用某种极限过程逐步逼近准确解的方法。
4.这次实践，将采用直接法：列主元高斯消去法，全主元高斯消去法；迭代法：雅克比迭代法和高斯-赛德尔迭代法。
三 详细设计 设有线性方程组11112211n n a x a x a x b +++=21122222n n a x a x a x b +++= 1122n n nn n n a x a x a x b +++=令A= 111212122212n n n n nn a a a aa a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ x=12n x x x ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ b= 12n b b b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦一. 上三角方程组的求解很简单。
所以若能将方程组化为上三角方程组，那就很容易得到方程的解，高斯消去法则是一种简洁高效的方法，可以将方程组化为上三角方程组，但是这种方法必须满足每一步时0kk a =才能求解，还有当kk a 绝对值很小时，作分母会引起较大的舍入误差。
因此在消元过程中应尽量选择绝对值较大的系数作为主元素。
1.列主元消去法很好的解决这一问题。
将1x 的系数1(1)k a k n ≤≤中绝对值最大者作为主元素，交换第一行和此元素所在的行，然后主元消去非主元项1x 的系数，。
实验五(线性方程组的数值解法和非线性方程求解)
1大学数学实验 实验报告 | 2014/4/5一、 实验目的1、学习用Matlab 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2、通过实例学习用线性代数方程组解决简化问题。
二、 实验内容项目一：种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变。
种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1,2,…,n ,当年年龄k 的种群数量记作x k ，繁殖率记作b k （每个雌性个体1年的繁殖的数量），自然存活率记作s k （s k =1−d k ，d k 为1年的死亡率），收获量记作ℎk ，则来年年龄k 的种群数量x ̌k 应该为x ̌k =∑b k n k=1x k , x ̌k+1=s k x k −ℎk , (k=1,2,…,n -1)。
要求各个年龄的种群数量每年维持不变就是要求使得x ̌k =x k , (k=1,2,…,n -1).(1) 如果b k , s k 已知，给定收获量ℎk ，建立求各个年龄的稳定种群数量x k 的模型（用矩阵、向量表示）.(2) 设n =5,b 1=b 2=b 5=0,b 3=5,b 4=3,s 1=s 4=0.4,s 2=s 3=0.6,如要求ℎ1~ℎ5为500,400,200,100,100，求x 1~x 5.(3) 要使ℎ1~ℎ5均为500，如何达到？问题分析：该问题属于简单的种群数量增长模型，在一定的条件（存活率，繁殖率等）下为使各年龄阶段的种群数量保持不变，各个年龄段的种群数量将会满足一定的要求，只要找到种群数量与各个参量之间的关系，建立起种群数量恒定的方程就可以求解出各年龄阶段的种群数量。
模型建立：根据题目中的信息，令x ̌k =x k ，得到方程组如下：{x ̌1=∑b k nk=1x k =x 1x ̌k+1=s k x k −ℎk =x k+1整理得到：{−x 1∑b k nk=1x k =0−x k+1+s k x k =ℎk2 大学数学实验 实验报告 | 2014/4/52写成系数矩阵的形式如下：A =[b 1−1b 2b 3s 1−100s 2−1…b n−1b n0000⋮⋱⋮000000000⋯00−10s n−1−1]令h =[0, ℎ1,ℎ2,ℎ3,…,ℎn−2,ℎn−1]Tx =[x n , x n−1,…,x 1]T则方程组化为矩阵形式：Ax =h ，即为所求模型。
数学实验练习整理(课本)
1. 统计推断（实验12）—区间估计、假设检验[mu,sigma,muci,sigmaci]=normfit(x,alpha); %%正态分布检验 [ht,sigt,cit]=ttest(x,mu); %%t 检验[hz,sigz,ciz,zval]=ztest(x,mu,sigma,alpha,tail); %%z 检验 tail 默认为0① P297第2题：（1）分别用两个月的数据验证这种说法的可靠性； 编程：x1=[]; x2=[]; alpha=0.05;[mu1,sigma1,muci1,sigmaci1]=normfit(x1,alpha) %%一月份的均值和标准差以及其置信区间 [mu2,sigma2,muci2,sigmaci2]=normfit(x2,alpha) %%二月份的均值和标准差以及其置信区间 运行结果：（1月）mu1 =115.1500； sigma1 =3.8699；muci1 =113.3388 116.9612； sigmaci1 = 2.9430 5.6523 （2月）mu2 =120.7500； sigma2 =3.7116muci2 =119.0129 122.4871； sigmaci2 =2.8227 5.4211（2）分别给出1月和2月汽油价格的置信区间（05.0=α）； 编程：x1=[]; x2=[]; mu=115; alpha=0.05;[h1,sigma1,ci1]=ttest(x1,mu,alpha,0) %%一月份汽油价格的置信区间 [h2,sigma2,ci2]=ttest(x2,mu,alpha,0) %%二月份汽油价格的置信区间 运行结果：（1月）h1 =0； sigma1 =0.8642； ci1 =113.3388 116.9612（2月）h2 =1； sigma2 =1.3241e-006； ci2 =119.0129 122.4871（3）如何给出1月和2月汽油价格差的置信区间（05.0=α） 编程：x1=[]; x2=[]; alpha=0.05;[h1,sigma1,ci1]=normfit(x2-x1,alpha) %数据看成同一个加油的数据，其价格差和置信区间 [h2,sigma2,ci2]=ttest(x2,x1,alpha,0) %数据完全随机时，用总体的t 分布检验 运行结果：h1 = 5.6000； sigma1 =5.4715； ci1 =3.0393 8.1607 h2 =1； sigma2 =2.0582e-004； ci2 =3.0393 8.1607结果分析：根据运行结果，我们可以知道数据完全随机时，用t 分布检验获得的结果更为合理准确。
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】【题目1】通过求解线性方程组A1x=b1和A2x=b2，理解条件数的意义和方程组的性态对解的影响。
其中A1是n阶范德蒙矩阵，即⎡1x0⎢1x1⎢A1=⎢⎢⎢⎣1xn-12x0x12 2xn-1n-1⎤ x0⎥ x1n-1⎥1,...,n-1 ，xk=1+0.1k，k=0,⎥ n-1⎥ xn-1⎥⎦A2是n阶希尔伯特矩阵，b1，b2分别是A1，A2的行和。
（1）编程构造A1（A2可直接用命令产生）和b1，b2；你能预先知道方程组A1x=和A2x=。
b2的解吗？令n=5，用左除命令求解（用预先知道的解可检验程序）b1（2）令n=5,7,9,…，计算A1，A2的条件数。
为观察它们是否病态，做以下试验：b1，b2不变，A1和A2的元素A1(n,n)，A2(n,n)分别加扰动ε后求解；A1和A2不变，b1，b2的分量b1(n)，分析A和b的微小扰动对解的影响。
b2(n)分别加扰动ε求解。
ε取10-1010，-8，10-6。
（3）经扰动得到的解记做x~，计算误差-x~x，与用条件数估计的误差相比较。
1.1构造A1，A2和b1，b2首先令n=5，构造出A1，A2和b1，b2。
首先运行以下程序，输出A1。
运行以下程序对A1,A2求行和：由于b1，b2分别是A1，A2的行和，所以可以预知x1=运行下列程序，用左除命令对b1，b2进行求解：得到以下结果： T。
x2=(1,1, ,1)1.2 计算条件数并观察是否为病态1.不加扰动，计算条件数。
运行以下程序：由此可知，A1，A2的条件数分别是3.574∗10, 4,766∗10。
2.b1，b2不变，A1(n,n)，A2(n,n)分别加扰动（1）n=5时设x11,x12,x13分别为A1添加扰动10−10，10−8,10−6后的解。
线性方程组的数值解法及其应用
线性方程组的数值解法及其应用一、问题描述现实中的问题大多数是连续的，例如工程中求解结构受力后的变形，空气动力学中计算机翼周围的流场，气象预报中计算大气的流动。
这些现象大多是用若干个微分方程描述。
用数值方法求解微分方程（组），不论是差分方法还是有限元方法，通常都是通过对微分方程（连续的问题，未知数的维数是无限的）进行离散，得到线性方程组（离散问题，因为未知数的维数是有限的）。
因此线性方程组的求解在科学与工程中的应用非常广泛。
经典的求解线性方程组的方法一般分为两类：直接法和迭代法。
二、基本要求1）掌握用MATLAB软件求线性方程初值问题数值解的方法；2）通过实例学习用线性方程组模型解决简化的实际问题；3）了解用高斯赛德尔列主元消去法和雅可比迭代法解线性方程组。
三、测试数据1） 直接法：A=[0.002 52.88;4.573 -7.290];b=[52.90;38.44];2） 迭代法：A=[10 -1 -2;-1 10 -2;-1 -1 5];b=[7.2;8.3;4.2];四、算法程序及结果1）function[RA,RB,n,x]=liezy1(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('因为RA~=RB,所以此方程组无解.')returnif RA==RBif RA==ndisp('因为RA=RB=n,所以此方程组有唯一解.')x=zeros(n,1);C=zeros(1,n+1);for p=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);x(n)=b(n)/A(n,n);for q=n-1:-1:1x(q)=(b(q)-sum(A(q,q+1:n)*x(q+1:n)))/A(q,q);endelsedisp('因为RA=RB<n,所以此方程组有无穷多解.')endend测试：A=[0.002 52.88;4.573 -7.290];>> b=[52.90;38.44];>> [RA,RB,n,x]=liezy1(A,b)因为RA=RB=n,所以此方程组有唯一解.RA =2RB =2n =2x =10.00001.00002）function Jacobi(A,b,x0,P,error,max1)[n n]=size(A);x=zeros(n,1);for k=1:max1for j=1;nx(j)=(b(j)-A(j,[1:j-1,j+1:n])*x0([1:j-1,j+1:n]))/A(j,j);endxerrx=norm(x-x0,P);x0=x;x1=A\b;if(errx<error)disp('迭代次数k,精确解x1和近似解x分别是：')kx1xreturnendendif(errx>=error)disp('请注意：Jacobi迭代次数已经超过最大迭代次数max1.') end测试：A=[10 -1 -2;-1 10 -2;-1 -1 5];>>b=[7.2;8.3;4.2];>>x0=[0;0;0];>>Jacobi(A,b,x0,inf,0.001,100)n =3x =0.7200迭代次数k,精确解x1和近似解x分别是：k =2x1 =1.10001.20001.3000x =0.7200五、应用举例1）营养学家配制一种具有1200卡，30g蛋白质及300mg维生素C的配餐。
线性代数方程组的数值解法
) a
a / a (k1) kj
( k 1) kk
a a (k1)
ij
( k 1) ik
( j) kj
( j k 1,, n 1) (i k 1,, n; j k 1,n 1)
第n步：得到：
1
a (1) 12
a (1) 13
a (1) 1n
1
a (2) 23
a (2) 2n
a (1) 1n 1
A b 行初等变换 I x
相应地，计算公式可表述为：
对 k 1,2,, n, 依次计算
aaik((jjkk
) )
a(k 1) kj
a(k 1) ij
/
a(k kk
1)
a a (k 1) (k )
ik
kj
(
j
k
1, k
2,, n 1)
(i 1,2,, k 1, k 1, k 2,, n; j k 1,, n 1)
后用第i行元数（i
2,
,
n)减去第一行对应元素的a
(0) i1
倍，(i 2,, n),这样，a1(01)位置变为1，其余各行的第一
个元素变为0。
1
(A(0)
,
b(0)
)
a
(0) 21
a
(0) n1
a(0) 12
/
a(0) 11
a(0) 22
a(0) 13
/
a(0) 11
a(0) 23
a(0) 1n
0
a (1) 13
a (2) 23
a (1) 1n 1
a (2) 2n 1
记为
a (2) 33
a
(2) 3n 1
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验6 线性代数方程组的数值解法[实验目的]1. 1. 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2. 2. 通过实例学习用线性代数方程组解决简化的实际问题。
[实验内容]5-5 输电网络：一种大型输电网络可简化为图5.5（见书）所示电路，其中R 1，R 2，…，R n 表示负载电阻，r 1，r 2，…，r n 表示线路内阻，I 1，I 2，…，I n 表示负载上的电流。
设电源电压为V 。
(1)列出求各负载电阻R 1，R 2，…，R n 的方程；(2)设I 1=I 2=…=I n =I ，r 1=r 2=…=r n =r ，在r=1，I=0.5，V=18，n=10的情况下求R 1，R 2，…，R n 及总电阻R 0。
[问题分析、模型建立及求解](1) 设电源负极为电势为0，电阻R 1上对应节点电压为V 1，对于任意节点，根据KCL 定律列出方程：111++----=k k k k k k k k r V V r V V R V而k kk R V I =，可得： 111111)(++++--++-=k k k k k k k k k k k k R r IR r I r I R r I Ik=2，3，……，n-1；k=1时2221211R r IR r I I +-=，为与上式形式一致，化为22212111111)(R r IR r I r I r V I +--=-k=m （12-≤≤n m ）时 111111)(++++--+--+=m m m n m m m m m m m m R r IR r I r I R r I Ik=n 时n n n n n n n R r IR r I I -=--11设以上方程组的矩阵形式为：b AR =则 []Tn R R R R 21=Tn I I I r V I b ⎥⎦⎤⎢⎣⎡-= 3211⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---------=------n n nn n n nn n n n n r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I A 1111114443333233322221222111000000(2)代入参数：5.021=====I I I I n ，121=====r r r r n ，V=18，n=10，⇒⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-----=5.05.0005.015.005.015.005.015.0005.01A ⇒[]Tb 5.05.05.05.17 -=在命令窗口输入MA TLAB 程序如下：clear all ;n=10; %由题目要求设定A11=sparse(1:n-1,1:n-1,-1,n,n); %定义A 的对角元素，除（n,n) A12=sparse(n,n,-0.5,n,n); %定义(n,n)A1=A11+A12; %对角元素A2=sparse(1:n-1,2:n,0.5,n,n); %输入A 的上次对角元素 A3=sparse(2:n,1:n-1,0.5,n,n); %输入A 的下次对角元素 A=A1+A2+A3;b1=0.5*ones(n,1); %b 的除第一项元素 b2=sparse(1,1,18,n,1); %b 的第一项元素 b=b1-b2; R=A\b输出结果如下：R =26.000017.0000 9.0000 2.0000 -4.0000 -9.0000 -13.0000 -16.0000 -18.0000 -19.0000所以各阻值为(R 1，R 2，…，R 10)=(26，17，9，2，-4，-9，-13，-16，-18，-19)总电阻R 0（即输入等效电阻）为00I V R =，又nII I nk k ==∑=10得到 )(6.35.010180Ω=⨯=R5-6 有5个反应器连接如图5.6（见书），各个Q 表示外部输入、输出及反应器间的流量（m 3/min ），各个c 表示外部输入及反应器内某物质的浓度(mg/m 3)。
假定反应器内的浓度是均匀的，利用质量守恒准则建立模型，求出各反应器内的浓度c 1~c 5，并讨论反应器j 外部输入改变１个单位(mg/min)所引起的反应器i 浓度的变化。
[问题分析、模型建立及求解]当反应容器中反应物浓度稳定时，输入物质质量与输出物质质量平衡，即输入物质质量等于输出物质质量。
由此分析，分别对1~5容器列质量平衡方程，以输入物质为“+”，输出物质为“－”⎪⎪⎪⎩⎪⎪⎪⎨⎧=+-+=+-+-=+-=++--=++-0)(0)(0)()(555542251155544443342240303334312232252423112010133111215c Q Q c Q c Q c Q c Q c Q c Q c Q c Q Q c Q c Q Q Q c Q c Q c Q c Q Q ，代入数据整理为⎪⎪⎪⎩⎪⎪⎪⎨⎧=++=+-+-=+=--=+-0430211816090335065215432322131c c c c c c c c c c c c c由矩阵表示：b Ac =[]Tc c c 51 =，[]Tb 00160050--=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡---=40013211810009100003300106A编写MATLAB 程序如下：A=[-6,0,1,0,0 3,-3,0,0,0 0,1,-9,0,0 0,1,8,-11,2 3,1,0,0,-4]; b=[-50,0,-160,0,0]';c=A\b %计算并输出c pause;dQc01=[1,0,0,0,0]';dQc02=[0,0,1,0,0]';c11=A\(b+dQc01); %01的输入减少1个单位 dc11=c-c2;c12=A\(b-dQc01); %01的输入增加1个单位 dc12=c-c1;c21=A\(b+dQc02); %03的输入减少1个单位 dc21=c-c4;c22=A\(b-dQc02) ; %03的输入增加1个单位 dc22=c-c3;[c,c11,dc11,c12,dc12,c21,dc21,c22,dc22]%外部输入改变引起的 %反应器浓度变化列表比较MATLAB 给出结果如下：c =11.5094 11.5094 19.0566 16.9983 11.5094所以改变前各反应器浓度为(c 1，c 2，c 3，c 4，c 5)=( 11.5094，11.5094，19.0566，16.9983，11.5094) 改变容器01或容器03的外部输入，得到各容器平衡时浓度及其增量如下表：[改进做法]上面这种做法显得有些麻烦，由b Ac =可得b A c 1-=，表明容器平衡浓度c 对外部输入b 是线性的，所以当b 增加1个单位（记作b ∆）时，c 的增量为b A c ∆=∆-1，若外部输入增加一个单位，Tb )0,0,0,0,1(=∆时，c ∆为1-A 的第一列， Tb )0,0,1,0,0(=∆时，c ∆为1-A 的第三列。
外部输入减少1个单位时，取其负值即可。
故改用矩阵求逆的方法来计算：A=[-6,0,1,0,0 3,-3,0,0,0 0,1,-9,0,0 0,1,8,-11,2 3,1,0,0,-4]; b=[-50,0,-160,0,0]';dx=inv(A) %求A 的逆矩阵输出结果为：dx =-0.1698 -0.0063 -0.0189 0 0 -0.1698 -0.3396 -0.0189 0 0 -0.0189 -0.0377 -0.1132 0 0 -0.0600 -0.0746 -0.0875 -0.0909 -0.0455 -0.1698 -0.0896 -0.0189 0 -0.2500红色部分显示为所求，此结果与前面的方法计算的一致，但工作量明显少了许多。
5-8 种群的繁殖与稳定收获：种群的数量因繁殖而增加，因自然死亡而减少，对于人工饲养的种群（比如家畜）而言，为了保证稳定的收获，各个年龄的种群数量应维持不变，种群因雌性个体的繁殖而改变，为方便起见以下种群数量均指其中的雌性。
种群年龄记作k=1，2，…，n ，当年年龄k 的种群数量记作x k ，繁殖率记作b k （每个雌性个体1年繁殖的数量），自然存活率记作s k （s k =1-d k ，d k 为1年的死亡率），收获量记作h k ，则来年年龄k 的种群数量k x ~应为)1,,2,1(~,~11,1-=-==+=∑n k h x s x x b x k k k k nk k k ，要求各个年龄的种群数量每年维持不变就是要使),,1(~n k x x k k ==。
(1)若已知b k ，s k ，给定收获量h k ，建立求个年龄的稳定种群数量x k 的模型（用矩阵、向量表示）。
(2)设n=5，b 1=b 2=b 5=0，b 3=5，b 4=3，s 1=s 4=0.4，s 2=s 3=0.6，如果要求h 1~h 5为500，400，200，200，100，求 (3)要使h 1~h 5均为500， 如何达到？[问题分析、模型建立及求解](1)要使各年龄种群数量每年维持不变即),,1(~n k x x k k ==，依题意得⎪⎩⎪⎨⎧-=-==+=∑)1,,2,1(111n k h x s x x b x k k k k nk kk用矩阵形式表示原方程组为：h Ax =[]Tn x x x x 21=， []Tn h h h h 0121-=nn n n b b b b s s s A ⨯-⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----=3211211100010001(2)代入题中数据⇒⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-----=0350114.016.016.014.0A ， []T h 0100200400500=编写MATLAB 程序如下：format bank ;A=[0.4,-1,0,0,0 0,0.6,-1,0,00,0,0.6,-1,0 0,0,0,0.4,-1 -1,0,5,3,0];h=[500,400,200,100,0]'; x=A\hMATLAB 输出结果如下：x =8481.01 2892.41 1335.44 601.27 140.51第5年龄段：x 5=140.5>100=h 5 ，说明收获量h 5可以达到100。
(3)要使h 1~h 5均为500，则h 变为：[]h Ax h ==0500500500500其他参数不变，程序变为：format bank ;A=[0.4,-1,0,0,0 0,0.6,-1,0,0 0,0,0.6,-1,0 0,0,0,0.4,-1 -1,0,5,3,0];h=[500,500,500,500,0]'; x=A\hMATLAB 输出结果如下：x =10981.01 3892.41 1835.44 601.27 -259.49从结果看出，x 5为-259.49，但种群数量不可能为负数，在本题所给条件下，无法使h 1~h 5均为500。