实验6 线性代数方程组的数值解法
线性代数方程组的数值解法讨论
线性代数方程组的数值解法讨论解线性方程组的方法，主要分为直接方法和迭代方法两种。
直接法是在没有舍入误差的假设下能在预定的运算次数内求得精确解。
而实际上，原始数据的误差和运算的舍入误差是不可以避免的，实际上获得的也是近似解。
迭代法是构造一定的递推格式，产生逼近精确解的序列。
对于高阶方程组，如一些偏微分方程数值求解中出现的方程组，采用直接法计算代价比较高，迭代法则简单又实用，因此比较受工程人员青睐。
小组成员本着工程应用，讨论将学习的理论知识转变为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 分解法等，另一类是迭代法（近似法），包括雅克比迭代法、高斯迭代法、超松弛迭代法等。
线性方程组的数值解法实验报告
实验报告——线性方程组的数值解法姓名：班级：学号：日期：一 实践目的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 的系数，。
数学实验——线性代数方程组的数值解
数学实验——线性代数⽅程组的数值解实验5 线性代数⽅程组的数值解法分1 黄浩 43⼀、实验⽬的1.学会⽤MATLAB软件数值求解线性代数⽅程组，对迭代法的收敛性和解的稳定性作初步分析；2.通过实例学习⽤线性代数⽅程组解决简化的实际问题。
⼆、For personal use only in study and research; not forcommercial use三、四、实验内容1.《数学实验》第⼆版（问题1）问题叙述：通过求解线性⽅程组，理解条件数的意义和⽅程组性态对解的影响，其中是n阶范德蒙矩阵，即是n阶希尔伯特矩阵，b1，b2分别是的⾏和。
（1）编程构造（可直接⽤命令产⽣）和b1，b2；你能预先知道⽅程组和的解吗？令n=5，⽤左除命令求解（⽤预先知道的解可验证程序）。
（2）令n=5,7,9，…，计算和的条件数。
为观察他们是否病态，做以下试验：b1，b2不变，和的元素，分别加扰动后求解；和不变，b1，b2的分量b1(n),b2(n)分别加扰动后求解。
分析A与b的微⼩扰动对解的影响。
取10^-10,10^-8,10^-6。
（3）经扰动得到的解记做，计算误差，与⽤条件数估计的误差相⽐较。
模型转换及实验过程：（1）⼩题.由b1,b2为,的⾏和，可知⽅程组和的精确解均为n ⾏全1的列向量。
在n=5的情况下，⽤matlab编程（程序见四.1），构造,和b1，b2，使⽤⾼斯消去法得到的解x1,x2及其相对误差e1,e2（使⽤excel计算⽽得）为：由上表可见，当n=5时，所得的解都接近真值，误差在10^-12的量级左右。
（2）⼩题分别取n=5,7,9,11,13,15，计算和的条件数c1和c2，（程序见四.2），结果如下：由上表可见，⼆者的条件数都⽐较⼤，可能是病态的。
为证实和是否为病态，先保持b不变，对做扰动，得到该情况下的⾼斯消元解，（程序见四.3），结果如下：(为使结果清晰简洁，在此仅列出n=5,9,13的情况，n=7,11,15略去)=10^-10时：=10^-8时：=10^-6时：由上表可见：a)对于希尔伯特阵，随着阶数的增加，微⼩扰动对解带来的影响越来越⼤，到了n=9时，已经有了6倍误差的解，到了n=13时，甚⾄出现了22倍误差的解元素；⽽随着的增加，解的偏差似乎也有增加的趋势，但仅凭上述表格⽆法具体判断（在下⼀⼩题中具体叙述）。
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】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后的解。
数值分析实验报告-清华大学--线性代数方程组的数值解法
数值分析实验报告-清华大学--线性代数方程组的数值解法(总15页)--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--线性代数方程组的数值解法实验1. 主元的选取与算法的稳定性问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。
主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。
实验内容：考虑线性方程组 n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求：（1）取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元，结果如何？（2）现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。
（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。
（4）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。
重复上述实验，观察记录并分析实验结果。
程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A);Ab=[A b];r=input('选主元方式(0:自动;1:手动)，r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元，请输入第',num2str(i),'列所选元素所处的行数：']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx运行结果(1)n=10，矩阵的条件数及自动选主元Cond(A,1) =×103Cond(A,2) = ×103Cond(A,inf) =×103程序自动选择主元（列主元）a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=0b.计算结果x=[1，1，1，1，1，1，1，1，1，1]T(2)n=10，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2…（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[, , , , , , , , , ]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k 步选择主元处于第k+1行） 最终计算得x=[1，1，1，1，1，1，1，1，1，1]T(3)n=20，手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：2…（实际选择时，第k 步选择主元处于第k 行） 最终计算得x=[，，，，，，，，，，，，，，，，，，，]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动)，r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元，请输入第1列所选元素所处的行数：2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元，请输入第2列所选元素所处的行数：3…（实际选择时，第k步选择主元处于第k+1行）最终计算得x=[1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1，1]T(4)A分别为幻方矩阵，Hilbert矩阵，pascal矩阵和随机矩阵简要分析计算(1)表明：对于同一矩阵，不同范数定义的条件数是不同的；Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
_第六章_线性方程组的数值解法迭代法
b 1
a 11
b2
f
a 22 bn
a nn
x(k1) B0x(k)f
--------(5)
第四节 解线性方程组的迭代法
令：
0 0 0
L
a 21
0
0 A的下三角部分矩阵
a n1 a n 2 0
0
U
0
a12 0
a1n a2n
A的上三角部分矩阵
第三节 向量范数和矩阵范数
（2）范数的另一个简单例子是二维欧氏空间的长度
0M x2 y2
欧氏范数也满足三个条件：
（勾股定理）
设x = (x1, x2) ① x 0 x >0 ② ax = a x a为常数 ③ x+ y ≤ x + y 前两个条件显然，第三个条件在几何上解释为三角形一边的长度不大于其它 两边长度之和。因此，称之三角不等式。
满足：
① A0,且A0,当且A 仅 0当
，若 A
正定
② A A,为任意实数
奇次
③ ABAB,A和 B为任意 n阶两 方个 三阵 角不等
则称 A 为矩阵A的范数。
第三节 向量范数和矩阵范数
2、矩阵范数与向量范数的相容性 对于任意的n维向量x，都有：
Ax A x
这一性质称为矩阵范数与向量范数的相容性。
n
A
max
1in
j1
aij
A的每行绝对值之和的最大值, 又称A的行范数
第三节 向量范数和矩阵范数
（3）矩阵的2范数
2范数 ||A|2 | ： (AT A )
(AAT) ?
矩阵的谱半径：
矩阵B的诸特征值为： i(i1,2, ,n)
线性代数方程组的数值解法
) 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
第三章 线性代数方程组的数值解
另外，矩阵求逆也是人们熟悉的求解线性方程 的一种方法。
[T ] [ A] [ B]
1
1
(3.5)
2 但求逆矩阵 [ A] 时，要计算 n个( n -1)阶行列式，
和一个 n 阶行列式，它的计算工作量也是相当可 观的，对于n 较大的情况也没有什么现实意义。
现在计算机上常用的直接解法大多数是以 系数矩阵的三角形化为基础的。就是说， 先对方程组进行变换，使其化为等价的 （即具有相同解的）三角形方程组。由于 三角形方程组的求解十分容易，原方程的 求解问题即告解决。
假定把要求的n 阶线性的方程组(3.1)改写成如 下形式：
a (1)T1 a (1)T2 a (1)T3 a (1)Tn b (1) 11 12 13 1n 1 (1) (1) (1) (1) (1) a21 T1 a22 T2 a23 T3 a2 T b n n 2 (3.10) (1) (1) (1) (1) (1) an1 T1 an 2T2 an3 T3 ann Tn bn
T1 b1 / l11 Ti (bi li1T1 li 2T2 li ,i 1 Ti 1 ) / lii i 2, 3, , n (i j )
这个计算过程通常也只作前推过程。
三角形方程组的求解是很简单的。方程组(3.6)的 计算公式可归结为：
(3.8)
21111??????????????????????nijkjabaabbnjkakijkkkkkkikkikijkkkjj110111313iiai系数矩阵中对角线上的元素都不应为零因为在消元的过程中不断用0adet用用gauss消去法解线性代数方程组时为了能求到最后的解尽管已经具备了的条件并使解尽可能的精确应消去法解线性代数方程组时为了能求到最后的解尽管已经具备了的条件并使解尽可能的精确应注意如下两点
- 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 I R 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 nn 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 ， []Th 0100200400500=编写MATLAB 程序如下：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,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 变为：[]hAx 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。