实验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消去法解线性代数方程组时为了能求到最后的解尽管已经具备了的条件并使解尽可能的精确应消去法解线性代数方程组时为了能求到最后的解尽管已经具备了的条件并使解尽可能的精确应注意如下两点
实验5_线性代数方程组的数值解法
实验5 线性代数方程组的数值解法化工系 分0班 2010011805 张亚清【实验目的】1、 学会用MATLAB 软件数值求解线性代数方程组，对迭代法的收敛性和解的稳定性作初步分析；2、 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】 1、 题目3已知方程组Ax=b ，其中2020⨯ℜ∈A ，定义为试通过迭代法求解此方程组，认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求：（1）选取不同的初始向量)0(x和不同的方程组右端向量b ，给出迭代误差要求，用雅克比迭代法和高斯-赛德尔迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛，记录迭代次数，分析计算结果并得出你的结论；（2）取定右端向量b 和初始向量)0(x，将A 的主对角线元素成倍增长若干次，非主对角线元素不变，每次用雅克比迭代法计算，要求迭代误差满足5)()1(10-∞+<-k k x x ，比较收敛速度，分析现象并得出你的结论。
【问题分析】对于线性方程组Ax=b ，满足0≠ii a （i=1,2,...,n ）,将A 分解为A=D-L-U ，则雅克比迭代公式等价于如下的矩阵形式,...2,1,0,)(1)(1)1(=++=--+k b D x D L D x k k或b D f U L D B f x B xJ J J k J k 11)()1(),(,--+=+=+=。
类似的，Ax=b 的高斯-赛德尔迭代公式等价于如下矩阵形式b L D f U L D B f x B x J S G S G k S G k 11)()1()(,)(,-----+-=-=+=。
【问题解答】（1）选取初始向量T x )1,...,1,1()0(=，T b )1,...,1,1(=，迭代要求为4)()1(10-∞+<-k k x x 。
将A 按A=D-L-U 分解为如下三个矩阵：①对方程组进行雅克比迭代，利用MATLAB 编程计算。
06实验六 方程及方程组的解
也可以通过初等变换方法，用函数 求解， 也可以通过初等变换方法，用函数reff求解，rref是将一个 求解 是将一个 矩阵化成行最简形。 矩阵化成行最简形。 C=[A,B]; %C为方程组的增广矩阵 为方程组的增广矩阵 R=rref(C) ↙%将C化成行最简形 将 化成行最简形 R= 1.0000 0 0 0 0 0 1.0000 0 0 0 0 0 1.0000 0 0 0 0 0 1.0000 0 0 0 0 0 1.0000 的最后一列元素就是所求之解。 则R的最后一列元素就是所求之解。 的最后一列元素就是所求之解
实验六方程及方程组的解实验61线性方程组的解实验62非线性方程的解一二分法二简单迭代法三牛顿迭代法四应用举例一求线性方程组的唯一解或特解二求线性齐次方程组的通解三求非齐次线性方程组的通解实验61线性方程组的解我们将线性方程组的求解分为两类
实验六 方程及方程组的解
实验6.1 线性方程组的解 实验
一、求线性方程组的唯一解或特解 二、求线性齐次方程组的通解 三、求非齐次线性方程组的通解
解：用矩阵除法求解如下： 用矩阵除法求解如下： A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; B=[1 4 0]'; X=A\B ↙ Warning: Rank deficient, rank = 2, tol = 8.8373e-015. X= 0 0 -0.5333 0.6000
5
2.2662 -1.7218 1.0571 -0.5940 0.3188
x1 + x2 − 3 x3 − x4 = 1 的一个特解。 例2 求方程组 3 x1 − x2 − 3 x3 + 4 x4 = 4 的一个特解。 x + 5x − 9x − 8x = 0 2 3 4 1
实验6线性代数方程组的数值解法
实验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ΛΛΛO O M O O O O M M O O O O O M M O M O ΛΛΛ(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.01ΛΛΛO M O O O O M M O O O O O MM O M O ΛΛΛA⇒[]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)。
计算方法线性方程组数值解法
d
2
a3b3c3
x3
d3
an
1bn1cn
1
xn
1
d
n
1
anbn xn dn
其系数矩阵为三对角形，元素满足以下条件：
|b1|>|c1|>0
|bi|≥|ai|+|ci|，且aici≠0 i=2，3，……n-1; |bn|≥|an|>0。
可以采用追赶法求解
4
线性代数方面的计算方法就是研究求解线 性方程组的一些数值解法与研究计算矩阵 的特征值及特征向量的数值方法。
5
设有线性方程组
a11x1 a12x2 a1nxn b1 a21x1a22x2a2nxnb2 an1x1 an2x2 annxn bn
式中，aij，bi为已知常数，xi为待求的未知量。记
u
2
2
u 2 n
u n 1,n 1u n 1,n
u n n
10
若uii≠0（i=1,2,……n）,则由下至上依次回代得
xn yn / unn
xn1 ( yn1 xi yi
un1,n xn ) / un1,n1
n
uij x j ) / uii
0
a
( 2
2 2
)
a
( 2
2) ,k 1
a
( 2
2) ,k
a
( 2
2) ,n
a
( 2
2) ,n 1
0 A(k)
0 0
a
( k
k) ,k
a
( k
k) ,k 1
a
k
k ,n
a
( k
k) 1,n
1
应用数值分析线性代数方程组数值解法
线性代数方程组数值解法——直接法一、n 阶线性代数方程组b x A =⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n n n b b b x x x a a a a a a a a a 2121212222111211二、上（下）方程组与回代（前推）过程 1、上三角方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n b b b x x x a a a a a a 2121222112112、回代过程⎪⎩⎪⎨⎧-=⎪⎭⎫ ⎝⎛-==∑+=)1,2,1(,1 n i a x a b x a b x ii ni j j ij ii nn n n 3、下三角方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n n nn n n b b b x x x a a a a a a 2121212221114、前推过程⎪⎩⎪⎨⎧=⎪⎭⎫ ⎝⎛-==∑-=),2,1(,111111n i a x a b x a b x ii i j j ij ii 三、顺序Gauss 消去法 1、消元过程（1）n 阶方程组经1-k 次消元后，得到的新的系数矩阵)(k A 具有如下形式：⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=)()()()()2(2)2(22)1(1)1(12)1(11)(k nn k nk k knk kk n n k a a a a a a a a a A（2）在第k 次，计算乘数),1(,)()(n k i a a m k kk k ik ik +==，则)1(+k A元素的计算公式为：⎩⎨⎧+=-=+=-=++),1(,),1,(,)()()1()()()1(n k i b m b b n k j i a m a a k k ik k i k ik kj ik k ij k ij （3）完成1-n 次消元后：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡)()2(2)1(121)()2(2)2(22)1(1112111n n n n nn n nb b b x x x a a a a a a ）（）（ 2、回代过程⎪⎩⎪⎨⎧-=⎪⎭⎫ ⎝⎛-==∑+=)1,2,1(,)(1)()()()( n k a x a b x a b x k kk n k j j k kj k k k n nn n n n3、计算量估计（1）第k 次消元需要进行)1)((+--k n k n 次乘法和)(k n -次除法：6523)1()()(231111nn n k n k n k n n k n k -+=+--+-∑∑-=-= （2）回代过程乘除法次数：22)(211nn n k n n k +=+-∑-= （3）乘数法总次数：3323n n n MD -+=（4）加减法总次数：652323n n n AS -+=4、要求：系数矩阵A 的所有顺序主子式),1,2,1(0n n k k -=≠∆ 四、列主元Gauss 消去法1、对于非奇异矩阵A ，在第一步消元时，即使A 的第一个元素为零，但是A 的第一列元素中至少有一个不为零，把这个不为零的元素所在的行与第一行交换位置，即可进行消元。
数学实验报告——利用MALTAB计算线性代数方程组的数值解法
实验三线性代数方程组的数值解法一、迭代法求解方程组㈠问题描述给定方程组的矩阵A，通过迭代法求解方程组。
1、选取不同的初始向量和不同的右端项向量，给定误差要求，用两种迭代法计算；2、去顶右端项向量和初始向量，将A的主对角线元素成倍增长若干次，非主对角线元素不变，用雅克比迭代法计算。
㈡方法与公式1、雅克比迭代法2、高斯-赛德尔迭代法㈢结果与分析1、不同初始向量、不同右端项向量、不同精度要求(1)初始向量定为zeros(n,1);①b=zeros(n,1)迭代次数为0，直接得到结果。
③b = 1:n(2)初始向量定为one s(n,1)①b=zeros(n,1)事实上，迭代100次时，所得结果约为10^-32，已经可以认为是0，但是由于没有达到精度要求，故不算收敛。
②b=ones(n,1)④b = n:1(3)初始向量定为1:n①b=zeros(n,1)②b=ones(n,1)(4)初始向量定为n:1①b=zeros(n,1)②b=ones(n,1)④b = n:1(5)简要小结a.在个别情况下雅可比迭代法收敛速度极慢，但事实上没有达到收敛时其计算结果已经可以接受；b.要求的精度越高，迭代的次数越多，迭代的次数与所要求的精度的对数值近似呈线性，也就是说两者近似呈指数关系；b.高斯-赛德尔迭代法有着比雅可比更好的迭代特性；2、更改A的主对角线元素(1) b =20:1;初值 1:20(2) b =20:1;初值 20:1(3) b =1:20;初值 20:1(4) b =[3;5;2;6;8;23;5;8;32;63;23;5;2;12;0;23;1;564;2;65]; 初值 ones(20,1)(5)简要小结a.迭代的次数随着对角线元素的成倍的增长而降低，趋于一稳定值；b.右端项以及迭代初值仅当对角线元素较小时对迭代次数起有作用，对角线元素成倍数增加后，迭代次数不变。
3、总结由以上各个对比可以得出以下结论：a.使用迭代法求解方程组时时，要求的精度越高，迭代次数越大；b.高斯-赛德尔迭代法的迭代次数要比雅可比迭代法迭代次数低；c.雅可比迭代的次数随着矩阵A对角线元素的成倍的增长而降低，d.当矩阵A的对角线元素足够大时，雅可比迭代法的迭代次数趋于稳定值；㈣程序清单1、第一问中的雅可比迭代function [y,k] = jacobi(A,b,m,tol)D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);y = ones(n,1);BJ=D\(L+U);fJ=D\b;k=0;while norm(A*y-b)/norm(b)>tol && k<mk = k+1;y = BJ*y+fJ;end2、高斯-赛德尔迭代法function [y,k] = GuassSeidel(A,b,m,tol)D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);BG = (D-L)\U;FG = (D-L)\b;y = zeros(n,1);k=0;while norm(A*y-b)/norm(b)>tol && k<mk = k+1;y = BG*y+FG;end3、第二问中的雅可比迭代function [y,k] = jacobi1(A,b,m,tol) D = diag(diag(A));L = - tril(A,-1);U = - triu(A,1);n = length(A);yk = ones(20,1);BJ=D\(L+U);fJ=D\b;k=0;yk1 = BJ*yk+fJ;k = k+1;ttol = norm(yk1-yk,inf);while ttol>tol && k<mk = k+1;yk = yk1;yk1 = BJ*yk+fJ;ttol = norm(yk1-yk,inf);endy=yk1;4、第一问脚本n=20;A1 = sparse(1:n,1:n,3,n,n);A2 = sparse(1:n-1,2:n,-1/2,n,n);A3 = sparse(2:n,1:n-1,-1/2,n,n);A4 = sparse(1:n-2,3:n,-1/4,n,n);A5 = sparse(3:n,1:n-2,-1/4,n,n);A = A1+A2+A3+A4+A5;b =zeros(20,1);m=10000;for i = 3:10tol = 10^(-i);[xJ,k1(i)] = jacobi(A,b,m,tol); [xG,k2(i)] = GuassSeidel(A,b,m,tol); endk1k25、第二问脚本n=20;k = 1;A1 = sparse(1:n,1:n,3,n,n);A2 = sparse(1:n-1,2:n,-1/2,n,n);A3 = sparse(2:n,1:n-1,-1/2,n,n);A4 = sparse(1:n-2,3:n,-1/4,n,n);A5 = sparse(3:n,1:n-2,-1/4,n,n);b =[3;5;2;6;8;23;5;8;32;63;23;5;2;12;0;23;1;564;2;65]; m=10000;tol = 10^(-5);for k = 1:10A = k*A1+A2+A3+A4+A5;[xJ,k1(k)] = jacobi1(A,b,m,tol);endk1二、一年生植物的繁殖㈠问题描述对于5.1.2的假设，给定参数，试分别用追赶法、稀疏系数矩阵和满矩阵求解；若b 有10%误差，估计对结果的影响。
计算方法实验报告-线性方程组的数值解法
重庆大学学生实验报告实验课程名称计算方法开课实验室DS1421学院年级专业学生姓名学号开课时间至学年第学期1．实验目的（1）高斯列主元消去法求解线性方程组的过程（2）熟悉用迭代法求解线性方程组的过程（3）设计出相应的算法，编制相应的函数子程序2．实验内容分别用高斯列主元消去法 ,Jacobi 迭代法,Gauss--Saidel 迭代法,超松弛迭代法求解线性方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡------725101391444321131243301024321x x x x 3.实验过程解：（1）高斯列主元消去法编制高斯列主元消去法的M 文件程序如下：%高斯列主元消元法求解线性方程组Ax=b%A 为输入矩阵系数,b 为方程组右端系数%方程组的解保存在x 变量中format long;%设置为长格式显示,显示15位小数A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]b=[10,5,-2,7]'[m,n]=size(A);%先检查系数正确性if m~=nerror('矩阵A 的行数和列数必须相同');return;endif m~=size(b)error('b 的大小必须和A 的行数或A 的列数相同');return;end%再检查方程是否存在唯一解if rank(A)~=rank([A,b])error('A 矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');return;endc=n+1;A(:,c)=b; %(增广)for k=1:n-1[r,m]=max(abs(A(k:n,k))); %选主元m=m+k-1; %修正操作行的值if(A(m,k)~=0)if(m~=k)A([k m],:)=A([m k],:); %换行endA(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去 endendx=zeros(length(b),1); %回代求解x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k);enddisp('X=');disp(x);format short;%设置为默认格式显示,显示5位运行，结果如下所示：（2）Jacobi迭代法编制迭代计算的M文件程序如下:%Jacobi迭代法求解% A为方程组的增广矩阵clc;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7] MAXTIME=50;%最多进行50次迭代eps=1e-5;%迭代误差[n,m]=size(A);x=zeros(n,1);%迭代初值y=zeros(n,1);k=0;%进入迭代计算disp('迭代过程X的值情况如下:')disp('X=');while 1disp(x');for i=1:1:ns=0.0;for j=1:1:nif j~=is=s+A(i,j)*x(j);endy(i)=(A(i,n+1)-s)/A(i,i);endendfor i=1:1:nmaxeps=max(0,abs(x(i)-y(i))); %检查是否满足迭代精度要求 endif maxeps<=eps%小于迭代精度退出迭代for i=1:1:nx(i)=y(i);%将结果赋给xendreturn;endfor i=1:1:n%若不满足迭代精度要求继续进行迭代x(i)=y(i);y(i)=0.0;endk=k+1;if k>MAXTIME%超过最大迭代次数退出error('超过最大迭代次数,退出');return;endend运行该程序结果如下:（3）Gauss--Saidel迭代法编制求解程序Gauss_Seidel.m如下:%Gauss_Seidel.m% A为方程组的增广矩阵clc;format long;A=[2,10,0,-3,10;-3,-4,-12,13,5;1,2,3,-4,-2;4,14,9,-13,7][n,m]=size(A);%最多进行50次迭代Maxtime=50;%控制误差Eps=10E-5;%初始迭代值x=zeros(1,n);disp('x=');%迭代次数小于最大迭代次数,进入迭代for k=1:Maxtimedisp(x);for i=1:ns=0.0;for j=1:nif i~=js=s+A(i,j)*x(j);%计算和endendx(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值end%因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件if sum((x-floor(x)).^2)<Epsbreak;end;end;X=x;disp('迭代结果:');Xformat short;运行结果如下所示：（4）超松弛迭代法编写函数M文件如下：%SOR法求解%w=1.45%方程组系数矩阵clc;A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]%方程组右端系数b=[10,5,-2,7]'w=1.45;%最大迭代次数Maxtime=100;%精度要求Eps=1E-5;%以15位小数显示format long;n=length(A);k=0;%初始迭代值x=ones(n,1);y=x;disp('迭代过程:');disp('x=');while 1y=x;disp(x');%计算过程for i=1:ns=b(i);for j=1:nif j~=is=s-A(i,j)*x(j);endendif abs(A(i,i))<1E-10 | k>=Maxtimeerror('已达最大迭代次数或矩阵系数近似为0,无法进行迭代'); return;ends=s/A(i,i);x(i)=(1-w)*x(i)+w*s;endif norm(y-x,inf)<Eps%达到精度要求退出计算break;endk=k+1;enddisp('最后迭代结果:');%最后的结果X=x'%设为默认显示格式format short; 结果如下：4．实验环境及实验文件存档名实验环境：Matlab7.0文件存档名：Gauss.m，Jacobi.m，Gauss_Seidel.m，SOR.m5.实验结果及分析=1.0000=2.0000=3.0000=4.0000经过验证，高斯列主元消结果正确。
- 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。