《大学数学实验》作业
线性方程组的数值解法
班级: 姓名: 学号: 日期:
目录
目录 .................................................................................................................................................. 2 【实验目的】 ................................................................................................................................... 3 【实验内容】 ................................................................................................................................... 3
【题目1】(课本习题第五章第3题) .................................................................................. 3
【第(1)问求解】 .............................................................................................................. 3 【第(2)问求解】 ............................................................................................................ 11 【本题小结】 ................................................................................................................. 12 【题目2】(课本习题第五章第5题) ................................................................................ 13
【第(1)问求解】 ............................................................................................................ 14 【第(2)问求解】 ............................................................................................................ 15 【本题小结】 ................................................................................................................. 16 【题目3】(课本习题第五章第7题) ................................................................................ 16
【第(1)问求解】 ....................................................................................................... 16 【第(2)问求解】 ....................................................................................................... 18 【本题小结】 ................................................................................................................. 19
【实验心得、体会】 ..................................................................................................................... 19
注:
本实验作业matlab脚本文件均以ex5_3形式命名,其中ex代表作业,5_3表示第五章第三小题
【实验目的】
1.学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;
2.通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】
【题目1】(课本习题第五章第3题)
已知方程组Ax=b,其中AR
2020
,定义为
1/21/43
1/231/21/4
1/41/231/21/4
1/41/23
A1/4
31/21/4
1/231/21/41/23
试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数
矩阵性质对收敛速度的影响。实验要求:
(1)选取不同初始向量x和不同的方程组右端向量b,给定迭代误差要求,用雅(0)
克比和高斯——赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;
(2)取定右端向量b和初始向量x(0),将A的主对角元素成倍增长若干次,非主
(k1)(k)5|xx||10对角元素不变,每次用雅克比迭代法,要求迭代误差满足,
比较收敛速度,分析现象并得出你的结论。
【第(1)问求解】
【模型分析】
对于矩阵A,将其分解为:
A=D-L-U,
其中:
30D=
00
0
000
30000
000, L=
030
0
0030
0
000
,U=
00
00
000
000
0
0
0
00 0
其中D、L、U均为20×20的方阵。
则雅可比迭代公式为
x(k1)D1(LU)x(k)D1b, k0,1,2...
高斯—赛德尔迭代公式为:
x(k1)D1(Lx(k1)Ux(k))D1b, k0,1,2...
选取初始向量x
(0)
=[0,0,00],b=[1,1,11],迭代误差为:x(k1)x(k)
105,
用雅可比迭代公法和高斯—赛德尔迭代法对方程组Ax=b进行迭代计算,在Matlab
中编写如下的程序代码: 【matlab程序】
%------------------------------作业题5_3脚本M文件源程序---------------------------- clear all;clc; n=20;
A1=sparse(1:n,1:n,3,n,n); % 输入A的对角元素 A2=sparse(1:n-1,2:n,-0.5,n,n); % 输入A的(上)次对角元素 A3=sparse(1:n-2,3:n,-0.25,n,n); % 输入(上)次次对角元素 A=A1+A2+A3+A2'+A3'; % 得到A矩阵 D=diag(diag(A)); % 提取D矩阵 L=-tril(A,-1); % 提取L矩阵 U=-triu(A,1); % 提取U矩阵
b= ones(20,1); % 设定b初值(进行变化) x(:,1)=zeros(n,1); % 设定初始向量(进行变化) for i=2:1000
x(:,i)=D\ (L+U)*x(:,i-1)+D\b; % 按雅克比迭代公式进行计算 if norm(x(:,i)-x(:,i-1),inf)
Bg_s= (D-L)\U; Fg_s= (D-L)\b;
y(:,1)=x(:,1); % 设定初始向量
for j=2:1000
y(:,j)=Bg_s*y(:,j-1)+Fg_s; % 按高斯—赛德尔迭代公式进行计算 if norm(y(:,j)-y(:,j-1),inf)
[i-1,j-1] % 输出两次迭代的次数
[x(:,i),y(:,j), A\b] % 输出迭代得到的方程的解及精确解 x,y % 输出两次迭代的所有结果
plot(1:i,x(1,:),'r',1:j,y(1,:),'b'); % 作出x1在两种迭代过程中的变化图像 xlabel('迭代次数');ylabel('x1');title('x1在两种迭代过程中的变化图像'); % 加入X轴、Y轴标记和标题
legend('雅克比迭代法','高斯—赛德尔迭代法', 'location', 'south') % 加入图例 gtext('x1');
可以得到,所得的向量序列均是收敛的。具体的原因见后面的理论分析部分。 用雅可比迭代公式计算,共迭代16次,得到的迭代结果如下表一:
表1:取x
(0)
=[0,0,00],b=[1,1,11]雅可比迭代公式计算部分结果
用高斯—赛德尔迭代法计算,共迭代11次,得到的部分迭代结果如下表2:
表2:取x
(0)
=[0,0,00],b=[1,1,11]高斯—赛德尔迭代法计算部分结果
x(0)=[0,0,00],b=[1,1,11],x(k1)x(k)
105时,为了更清晰比较两种迭代
方法收敛快慢,作出解x1随迭代次数变化的图像:
x1在两种迭代过程中的变化图像
x1
迭代次数
图1: 解x1随迭代次数变化图像
【结果分析】
由以上数据图表可知,x
(0)
=[0,0,00],b=[1,1,11],x(k1)x(k)
105时,解
向量序列均是收敛的,都收敛到精确解,不过收敛速度不同。雅克比迭代法共迭代16次,高斯—赛德尔迭代法共迭代11次,因此可以得出高斯—赛德尔迭代法收敛速度更快。两种迭代法的收敛情况类似。迭代前几步变化较大,后面均在小范围内变化,与真实解已十分接近,直到满足相应的迭代误差条件时停止迭代过程。这从图1可以很直观地看出。
选取不同初始向量x和不同的方程组右端向量b,给定不同迭代误差要求,比较雅克比和高斯——赛德尔迭代法
① 选取相同初始向量x和不同的方程组右端向量b,给定相同迭代误差要求,用雅克比和高斯——赛德尔迭代法计算,比较结果:
表3:
、迭代误差相同,不同右端向量b时两种迭代方法的计算比较
(0)(0)
【结果分析】
由表格可知,仅右端向量b不同时,解向量序列均是收敛的,在误差允许范围内,
对相同b,两种方法得到解向量相同;不同b之间,解向量不同。对相同b,两种方法收敛速度不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多,收敛速
度慢,高斯—赛德尔迭代法收敛速度快。还可以看出,b向量越复杂,值越大,两种迭代方法的迭代次数均增加。 分析:具体的原因见后面的理论分析部分。
② 选取不同初始向量x给定相同的方程组右端向量b,迭代误差要求,用雅克比和高斯——赛德尔迭代法计算,比较结果:
表4:
不同,向量b、迭代误差相同时两种迭代方法的计算比较
(0)
【结果分析】
由表格可知,仅初始向量x不同时,解向量序列均是收敛的,在误差允许范围内解向量均相同,不过收敛速度不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多,收敛速度慢,高斯—赛德尔迭代法收敛速度快。还可知,初始向量x越大,偏离解越多,所需要的迭代次数越多,两种迭代方法的迭代次数均增加。 分析:这从直观上也是很容易理解的:当初始
距离解x越远时,所需要的修
(0)
(0)
正次数就越多,故迭代次数越多。因为对于同一个b,方程形式是相同的,故在收敛的前提下,不同的
最终将趋于同一个值,也就是方程的真实解。
③ 选取相同初始向量x、方程组右端向量b,给定不同迭代误差要求,用雅克比和高斯——赛德尔迭代法计算,比较结果:
表5:
、向量b相同、迭代误差不同时两种迭代方法的计算比较
(0)
【结果分析】
由表格可知,仅迭代误差要求不同时,解向量序列均是收敛的,即便给定的条件
相同,两种方法得到的解向量也略有不同,如表格中红色区域所示,说明迭代精度要求不高时,解略有偏差。对于相同条件,两种迭代法收敛速度也不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多,收敛速度慢,高斯—赛德尔迭代法收敛速度快。还可以看出,随着迭代误差要求提高,两种迭代方法的迭代次数均增加。
分析:误差条件变严格后,原来的迭代次数所对应的解已无法满足新的误差条件,故需要更多的迭代次数以使解满足当前的误差条件。
【理论分析】
①收敛性的分析
两种迭代方法在求解线性方程组时,都能写成下列形式:
从而有
进而得到下列结果:
(1)对于雅可比迭代法法,B=D1(LU)。在Matlab中运行eig(D_inv*(L+U)语句,得到B=D1(LU)的全部特征值为-0.4893,-0.4579,-0.4079,-0.3425,-0.2658,-0.1825,-0.0978,-0.0162,0.0577,0.1205,0.1677,0.1700,0.1819,0.1914,0.2088,0.2106,0.2298,0.2317,0.2445,0.2454。其中绝对值的最大的值为:0.4893<1,即B的谱半径小于1,故对方程Ax=b进行雅克比迭代是收敛的。
(2)对于高斯—赛德尔迭代法,
。同样可以求得B的全部特征值
绝对值最大的值为0.1210<1,即B的谱半径小于1,故对方程Ax=b进行高斯—赛德尔迭代是收敛的。
可见,上述分析与之前的运行结果是相一致的。 ②收敛速度分析
由于矩阵的谱半径不超过它的任一种范数,所以若敛,且有误差估计式
,则迭代公式收
可见,
越小,序列
收敛越快。
对于雅可比迭代法,B=D1(LU)。利用Matlab的norm命令可以求得B的∞-范数为0.5000。对于高斯—赛德尔迭代法,
。同样可以求得B的
∞-范数为0.2802。即高斯—赛德尔迭代法对应的序列
值更小,故与之对应的
收敛更快,从而对于一定的迭代误差条件,高斯—赛德尔迭代法比雅可
比迭代法法的迭代次数要少。
可见,上述分析与之前的运行结果是相一致的。
【拓展对比】
本题采用的是迭代法求解线性方程组。也可以通过直接法求解线性方程组,
如高斯消元法,LU分解等。可以在Matlab中利用x=A\b或者利用递推关系式进行求解。求出来的结果与用迭代法求出来的解略有出入。这是由求解过程中的误差所致。
【第(2)问求解】
取定方程组右端向量b=[1,2,320],同时取定初始向量x(0)=[0,0,00],并将迭代误差设为为:x(k1)x(k)
105,将A的主对角线元素分别增长为6,
12,24,48,非主对角元素不变,每次用雅可比迭代法进行迭代计算,观察收敛速度的变化,得到结果如下面的表6:
表6: A主对角元素加倍时收敛速度的变化(b=
,x(0)=
【结果分析】
由表6可以明显看到,随着矩阵A对角元素的增大,对于雅可比迭代法,与之对应的
的值越来越小。在b和
相同以及迭代误差条件相同的前提
下,迭代次数越来越少,即收敛速度越来越快。
分析: 由前一问对收敛速度的分析可得, 由于矩阵的谱半径不超过它的任一种范数,所以若收敛,且有误差估计式
,则迭代公式
可见,
越小,序列
收敛越快。
对于雅可比迭代法,B=D1(LU)。利用Matlab的norm命令可以求得当 A的主对角元素分别为3,6,12,24,48时,相应的B的∞-范数分别为为0.5000,0.2500,0.1250,0.0625,0.03125。即随着A的主对角元素越来越大,相应的的值越来越小,从而序列
收敛越快,故迭代次数越来越少(在相同的b和
,
以及相同的误差条件的前提下)。
【本题小结】
(1)对于本题中的方程组Ax=b,选取不同的初始向量和不同的方程组右端向量b,所得到的迭代向量序列均是收敛的。这是由两种迭代法所对应的B的
谱半径均小于1所决定的。
(2)对于固定的和b,以及相同的迭代误差条件,两种迭代方法的迭代次数不同。高斯—赛德尔迭代法对应的迭代次数要比雅可比迭代法对应的迭代次数少。这是因为高斯—赛德尔迭代法对应的收敛更快。
(3)对于同一个b,不同的对应的迭代次数不同,且距离解x越远,所需的迭代次数越多。但最后得到的解是一样的。而对于不同的b,最后得到的解一般是不同的。
(3)对于同一个,相同b,给定不同的迭代误差要求,最后得到的解可能略有偏差。误差精度越高,所需的迭代次数越多。误差条件变严格后需要更多的迭代次数以使解满足当前的误差条件。
(5)通过数据表和图像还可以观察到,两种迭代法的收敛情况类似。迭代前几步变化较大,后面均在小范围内变化,与真实解已十分接近,直到满足相应的迭代误差条件时停止迭代过程。
(6)随着矩阵A对角元素的增大,在b和相同以及迭代误差条件相同的前提下,迭代次数越来越少,即收敛速度越来越快。这是由与之对应的的值越来越小所决定的。
值更小,
故与之对应的序列
【题目2】(课本习题第五章第5题)
设国民经济由农业、制造业和服务业三个部门构成,已知某年它们之间的投入产
出关系、外部需求、初始投入等如表7所示。
表7
(1) 如果今年对农业、制造业和服务业的外部需求分别为50、150、100亿元,
问这部门的总产出分别应为多少?
(2) 如果三个部门的外部需求分别增加一个单位,问它们的总产出分别增加多
少?
【第(1)问求解】
由表7得到三个部门的投入产出表如表8所示:
表8
xi
设有n个部门,记一定时期内第i个部门的总产出为投入为
xij
,其中对第j个部门的
,外部需求为
di
,则
xixijdi,i1,2,...,n.
j1n
根据投入系数
aij
的定义,有
xijaijxj,i,j1,2,...,n,
即
aij
是第j个部门的单位产出所需的第i个部门的投入。由此可得
xiaijxjdi,i1,2,...,n.
j1n
,需求向量
记投入系数矩阵
A(aij)nn
,产出向量
x(x1,x2,...,xn)T
d(d1,d2,...,dn)T
,则上式可写作xAxd,或
(IA)xd.
已知今年对农业、制造业和服务业的外部需求分别为50亿元,150亿元,100亿元,求这三个部门的总产出: 【matlab程序】
%------------------------------作业题5_5脚本M文件源程序---------------------------- clear all;clc;
A=[0.15 0.1 0.2;0.3 0.05 0.3;0.2 0.3 0]; % 设定系数矩阵A d=[50 150 100]'; % 给定外部需求 b=eye(3)-A;
x=b\d % 求解方程
dx=inv(b)
得到结果如下:
x 139.2801, 267.6056, 208.1377
T
1.34590.25040.3443
dx0.56341.26760.4930
0.43820.43041.2167
即三个部门的总产出应分别为:139.2801、267.6056、208.1377(亿元)。
【第(2)问求解】
当三个部门的外部需求分别增加1个单位时,求它们总产出相应增加的量: 由(IA)xd可得
x(IA)1d.
表明总产出x对外部需求d是线性的,所以当d增加1个单位(记作d)时,x
1
x(IA)d. 的增量为
于是从dx=inv(b)矩阵中便可得三个部门的外部需求分别增加一个单位时,各部
门总产出分别增加量。
1.34590.25040.3443
dx0.56341.26760.4930
0.43820.43041.2167
结果如下:
T1
d(1,0,0)(IA)x若农业的外部需求增加1个单位,,为的第1列
x=1.3459,0.5634,0.4382;即当农业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加1.3459, 0.5634, 0.4382单位。
若制造业的外部需求增加1个单位,d(0,1,0),x为(IA)的第2列
T
1
T
x=0.2504,1.2676,0.4304;即当制造业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加0.2504,1.2676,0.4304单位。
若服务业的外部需求增加1个单位,d(0,0,1),x为(IA)的第3列
T
1
T
x=0.3443,0.4930,1.2167;即当服务业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加0.3443,0.4930,1.2167单位。
T
【本题小结】
本题是对实际问题的数学建模,模型较简单,问题最后转化为线性方程组求解问题。因为方程并不复杂,我采用直接法求解。很快得到了答案。
【题目3】(课本习题第五章第7题)
输电网络:一种大型输电网络可简化为图5.5所示电路,其中R1,R2,…,Rn表示负载电阻,r1,r2,…,rn表示线路内阻,I1,I2,…,In表示负载上的电流。设电源电压为V。
(1)列出求各负载电流I1,I2,…,In的方程;
(2)设R1=R2=…=Rn=R,r1=r2=…=rn=r,在r=1,R=6,V=18,V=18,n=10的情况下求I1,I2,…,In及总电流I0。
图5.5
【第(1)问求解】
问题分析及模型建立:
IIiIi1Inr
通过对以上电路的分析可得,流过电阻i的电流ri。对于每一个
由Ri,ri+1,Ri1组成的闭合回路,根据基尔霍夫回路电压法,有
-IiRi+Iri+1ri+1Ii1Ri10。对于由电源V,R1,r1组成的回路,有
-V+I1R1Ir1r10
。
因此,可列方程如下:
I1R1(I1+I2+In)r1V0
-IRIR(III)r0112223n2
-I2R2I3R3(I3I4In)r30 -In-1Rn1InRnInrn0
整理后可得,
R11(1+)I+I+IV12nr1r1
R2R1
-I(1)I2I3In=0r1
r22
RR-2I(13)III=0
34n
r32r3
-Rn1I(1Rn)I0
n-1n
rn
rn变成矩阵的形式,即AXb,其中,
1+R1 1 1 … … 1 r1
R1R2r 12r 1 … 12 A= R
2r 3
… Rn1r 1RnnrnIV1
I
2rX
1I03
b0
,。 In-10In
最后列出求各负载电流I1,I2,…,In的方程如下:
R1I1+ 1 1 … … 1 r1V
1
r1
R2R1
r 1r 1 … 1 I20
22
R2
I30 r3
… In-10
Rn1Rn
1 In0 rnrn
【第(2)问求解】
由题意可得, R1=R2=…=Rn=R,r1=r2=…=rn=r,r=1,R=6,V=18, n=10因此,
矩阵A变作
7 1 1 … 1
06 7 1 … 1
6 矩阵b变作 b=0 A=
0
6 7
7 1 1 … 1 I1186 7 1 … 1 I02
整理方程得到 6 I3=0
In-10 6 7In0
18
其中n=10。
【Matlab程序】
%------------------------------作业题5_7脚本M文件源程序---------------------------- clear all;clc; n=10;
A1=sparse(1:n,1:n,7,n,n); A2=sparse(2:n,1:n-1,-6,n,n); A3=triu(ones(10,10),1);
A=A1+A2+A3; % 设定系数矩阵A b=[18,zeros(1,9)]'; % 设定列向量b x=[A\b]' % 求解各支路电流值 I0=sum(x) % 求解电流I0
得到如下结果:
表9:r=1,R=6,V=18,V=18,n=10时I1,I2,…,In,I0的值
所以要求的I1,I2,…,In及总电流I0值如表9所示。
【本题小结】
本题是电路问题的求解。可以采用差分方程求解,也可以用线性方程组求解。我采用线性方程组求解得到答案。因为问题中矩阵A并不复杂,因此采用直接法求解方程组,模型建好后很快便得到了答案。
【实验心得、体会】
本次实验作为数学实验课程的第三次作业,让我们掌握了线性方程组的数值解法相关基础知识,同时我进一步了解了matlab。本次实验第一个题目难度较大其余两题目都比较轻松。通过第一题,我深刻理解了雅克比和高斯——赛德尔迭代法,增长了许多知识。实验中最令我困惑的是为什么很简单的一个矩阵求解(线性方程组求解)问题要用雅克比和高斯——赛德尔迭代法求解,比如第一问中已知矩阵A与b,直接用x=A\b便可得解,我们却要绕一个大弯得到x的解,后来我逐渐明白了,直接法适用于求解中小型(n
我尝试着对每个问题用学到的不同的方法去解决(如对线性方程组分别用直接法中的高斯消元法和LU分解以进行求解,并和两种迭代方法进行比较),同时尝试着从不同的方面和角度去分析,并尽可能地发掘出更多值得去思考和分析的东西,在这个过程中不仅巩固了此次所学到的了线性方程组的数值解等相关知识,也不仅进一步熟练了Matlab的操作,更重要的是通过这次作业进行了一次思维上的极大锻炼,进一步熟悉了通过数学建模思想和数学软件解决一些和现实密切相关的问题。
总之,本次实验收获挺大的。
《大学数学实验》作业
线性方程组的数值解法
班级: 姓名: 学号: 日期:
目录
目录 .................................................................................................................................................. 2 【实验目的】 ................................................................................................................................... 3 【实验内容】 ................................................................................................................................... 3
【题目1】(课本习题第五章第3题) .................................................................................. 3
【第(1)问求解】 .............................................................................................................. 3 【第(2)问求解】 ............................................................................................................ 11 【本题小结】 ................................................................................................................. 12 【题目2】(课本习题第五章第5题) ................................................................................ 13
【第(1)问求解】 ............................................................................................................ 14 【第(2)问求解】 ............................................................................................................ 15 【本题小结】 ................................................................................................................. 16 【题目3】(课本习题第五章第7题) ................................................................................ 16
【第(1)问求解】 ....................................................................................................... 16 【第(2)问求解】 ....................................................................................................... 18 【本题小结】 ................................................................................................................. 19
【实验心得、体会】 ..................................................................................................................... 19
注:
本实验作业matlab脚本文件均以ex5_3形式命名,其中ex代表作业,5_3表示第五章第三小题
【实验目的】
1.学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;
2.通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】
【题目1】(课本习题第五章第3题)
已知方程组Ax=b,其中AR
2020
,定义为
1/21/43
1/231/21/4
1/41/231/21/4
1/41/23
A1/4
31/21/4
1/231/21/41/23
试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数
矩阵性质对收敛速度的影响。实验要求:
(1)选取不同初始向量x和不同的方程组右端向量b,给定迭代误差要求,用雅(0)
克比和高斯——赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;
(2)取定右端向量b和初始向量x(0),将A的主对角元素成倍增长若干次,非主
(k1)(k)5|xx||10对角元素不变,每次用雅克比迭代法,要求迭代误差满足,
比较收敛速度,分析现象并得出你的结论。
【第(1)问求解】
【模型分析】
对于矩阵A,将其分解为:
A=D-L-U,
其中:
30D=
00
0
000
30000
000, L=
030
0
0030
0
000
,U=
00
00
000
000
0
0
0
00 0
其中D、L、U均为20×20的方阵。
则雅可比迭代公式为
x(k1)D1(LU)x(k)D1b, k0,1,2...
高斯—赛德尔迭代公式为:
x(k1)D1(Lx(k1)Ux(k))D1b, k0,1,2...
选取初始向量x
(0)
=[0,0,00],b=[1,1,11],迭代误差为:x(k1)x(k)
105,
用雅可比迭代公法和高斯—赛德尔迭代法对方程组Ax=b进行迭代计算,在Matlab
中编写如下的程序代码: 【matlab程序】
%------------------------------作业题5_3脚本M文件源程序---------------------------- clear all;clc; n=20;
A1=sparse(1:n,1:n,3,n,n); % 输入A的对角元素 A2=sparse(1:n-1,2:n,-0.5,n,n); % 输入A的(上)次对角元素 A3=sparse(1:n-2,3:n,-0.25,n,n); % 输入(上)次次对角元素 A=A1+A2+A3+A2'+A3'; % 得到A矩阵 D=diag(diag(A)); % 提取D矩阵 L=-tril(A,-1); % 提取L矩阵 U=-triu(A,1); % 提取U矩阵
b= ones(20,1); % 设定b初值(进行变化) x(:,1)=zeros(n,1); % 设定初始向量(进行变化) for i=2:1000
x(:,i)=D\ (L+U)*x(:,i-1)+D\b; % 按雅克比迭代公式进行计算 if norm(x(:,i)-x(:,i-1),inf)
Bg_s= (D-L)\U; Fg_s= (D-L)\b;
y(:,1)=x(:,1); % 设定初始向量
for j=2:1000
y(:,j)=Bg_s*y(:,j-1)+Fg_s; % 按高斯—赛德尔迭代公式进行计算 if norm(y(:,j)-y(:,j-1),inf)
[i-1,j-1] % 输出两次迭代的次数
[x(:,i),y(:,j), A\b] % 输出迭代得到的方程的解及精确解 x,y % 输出两次迭代的所有结果
plot(1:i,x(1,:),'r',1:j,y(1,:),'b'); % 作出x1在两种迭代过程中的变化图像 xlabel('迭代次数');ylabel('x1');title('x1在两种迭代过程中的变化图像'); % 加入X轴、Y轴标记和标题
legend('雅克比迭代法','高斯—赛德尔迭代法', 'location', 'south') % 加入图例 gtext('x1');
可以得到,所得的向量序列均是收敛的。具体的原因见后面的理论分析部分。 用雅可比迭代公式计算,共迭代16次,得到的迭代结果如下表一:
表1:取x
(0)
=[0,0,00],b=[1,1,11]雅可比迭代公式计算部分结果
用高斯—赛德尔迭代法计算,共迭代11次,得到的部分迭代结果如下表2:
表2:取x
(0)
=[0,0,00],b=[1,1,11]高斯—赛德尔迭代法计算部分结果
x(0)=[0,0,00],b=[1,1,11],x(k1)x(k)
105时,为了更清晰比较两种迭代
方法收敛快慢,作出解x1随迭代次数变化的图像:
x1在两种迭代过程中的变化图像
x1
迭代次数
图1: 解x1随迭代次数变化图像
【结果分析】
由以上数据图表可知,x
(0)
=[0,0,00],b=[1,1,11],x(k1)x(k)
105时,解
向量序列均是收敛的,都收敛到精确解,不过收敛速度不同。雅克比迭代法共迭代16次,高斯—赛德尔迭代法共迭代11次,因此可以得出高斯—赛德尔迭代法收敛速度更快。两种迭代法的收敛情况类似。迭代前几步变化较大,后面均在小范围内变化,与真实解已十分接近,直到满足相应的迭代误差条件时停止迭代过程。这从图1可以很直观地看出。
选取不同初始向量x和不同的方程组右端向量b,给定不同迭代误差要求,比较雅克比和高斯——赛德尔迭代法
① 选取相同初始向量x和不同的方程组右端向量b,给定相同迭代误差要求,用雅克比和高斯——赛德尔迭代法计算,比较结果:
表3:
、迭代误差相同,不同右端向量b时两种迭代方法的计算比较
(0)(0)
【结果分析】
由表格可知,仅右端向量b不同时,解向量序列均是收敛的,在误差允许范围内,
对相同b,两种方法得到解向量相同;不同b之间,解向量不同。对相同b,两种方法收敛速度不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多,收敛速
度慢,高斯—赛德尔迭代法收敛速度快。还可以看出,b向量越复杂,值越大,两种迭代方法的迭代次数均增加。 分析:具体的原因见后面的理论分析部分。
② 选取不同初始向量x给定相同的方程组右端向量b,迭代误差要求,用雅克比和高斯——赛德尔迭代法计算,比较结果:
表4:
不同,向量b、迭代误差相同时两种迭代方法的计算比较
(0)
【结果分析】
由表格可知,仅初始向量x不同时,解向量序列均是收敛的,在误差允许范围内解向量均相同,不过收敛速度不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多,收敛速度慢,高斯—赛德尔迭代法收敛速度快。还可知,初始向量x越大,偏离解越多,所需要的迭代次数越多,两种迭代方法的迭代次数均增加。 分析:这从直观上也是很容易理解的:当初始
距离解x越远时,所需要的修
(0)
(0)
正次数就越多,故迭代次数越多。因为对于同一个b,方程形式是相同的,故在收敛的前提下,不同的
最终将趋于同一个值,也就是方程的真实解。
③ 选取相同初始向量x、方程组右端向量b,给定不同迭代误差要求,用雅克比和高斯——赛德尔迭代法计算,比较结果:
表5:
、向量b相同、迭代误差不同时两种迭代方法的计算比较
(0)
【结果分析】
由表格可知,仅迭代误差要求不同时,解向量序列均是收敛的,即便给定的条件
相同,两种方法得到的解向量也略有不同,如表格中红色区域所示,说明迭代精度要求不高时,解略有偏差。对于相同条件,两种迭代法收敛速度也不同。雅克比迭代法比高斯—赛德尔迭代法迭代次数多,收敛速度慢,高斯—赛德尔迭代法收敛速度快。还可以看出,随着迭代误差要求提高,两种迭代方法的迭代次数均增加。
分析:误差条件变严格后,原来的迭代次数所对应的解已无法满足新的误差条件,故需要更多的迭代次数以使解满足当前的误差条件。
【理论分析】
①收敛性的分析
两种迭代方法在求解线性方程组时,都能写成下列形式:
从而有
进而得到下列结果:
(1)对于雅可比迭代法法,B=D1(LU)。在Matlab中运行eig(D_inv*(L+U)语句,得到B=D1(LU)的全部特征值为-0.4893,-0.4579,-0.4079,-0.3425,-0.2658,-0.1825,-0.0978,-0.0162,0.0577,0.1205,0.1677,0.1700,0.1819,0.1914,0.2088,0.2106,0.2298,0.2317,0.2445,0.2454。其中绝对值的最大的值为:0.4893<1,即B的谱半径小于1,故对方程Ax=b进行雅克比迭代是收敛的。
(2)对于高斯—赛德尔迭代法,
。同样可以求得B的全部特征值
绝对值最大的值为0.1210<1,即B的谱半径小于1,故对方程Ax=b进行高斯—赛德尔迭代是收敛的。
可见,上述分析与之前的运行结果是相一致的。 ②收敛速度分析
由于矩阵的谱半径不超过它的任一种范数,所以若敛,且有误差估计式
,则迭代公式收
可见,
越小,序列
收敛越快。
对于雅可比迭代法,B=D1(LU)。利用Matlab的norm命令可以求得B的∞-范数为0.5000。对于高斯—赛德尔迭代法,
。同样可以求得B的
∞-范数为0.2802。即高斯—赛德尔迭代法对应的序列
值更小,故与之对应的
收敛更快,从而对于一定的迭代误差条件,高斯—赛德尔迭代法比雅可
比迭代法法的迭代次数要少。
可见,上述分析与之前的运行结果是相一致的。
【拓展对比】
本题采用的是迭代法求解线性方程组。也可以通过直接法求解线性方程组,
如高斯消元法,LU分解等。可以在Matlab中利用x=A\b或者利用递推关系式进行求解。求出来的结果与用迭代法求出来的解略有出入。这是由求解过程中的误差所致。
【第(2)问求解】
取定方程组右端向量b=[1,2,320],同时取定初始向量x(0)=[0,0,00],并将迭代误差设为为:x(k1)x(k)
105,将A的主对角线元素分别增长为6,
12,24,48,非主对角元素不变,每次用雅可比迭代法进行迭代计算,观察收敛速度的变化,得到结果如下面的表6:
表6: A主对角元素加倍时收敛速度的变化(b=
,x(0)=
【结果分析】
由表6可以明显看到,随着矩阵A对角元素的增大,对于雅可比迭代法,与之对应的
的值越来越小。在b和
相同以及迭代误差条件相同的前提
下,迭代次数越来越少,即收敛速度越来越快。
分析: 由前一问对收敛速度的分析可得, 由于矩阵的谱半径不超过它的任一种范数,所以若收敛,且有误差估计式
,则迭代公式
可见,
越小,序列
收敛越快。
对于雅可比迭代法,B=D1(LU)。利用Matlab的norm命令可以求得当 A的主对角元素分别为3,6,12,24,48时,相应的B的∞-范数分别为为0.5000,0.2500,0.1250,0.0625,0.03125。即随着A的主对角元素越来越大,相应的的值越来越小,从而序列
收敛越快,故迭代次数越来越少(在相同的b和
,
以及相同的误差条件的前提下)。
【本题小结】
(1)对于本题中的方程组Ax=b,选取不同的初始向量和不同的方程组右端向量b,所得到的迭代向量序列均是收敛的。这是由两种迭代法所对应的B的
谱半径均小于1所决定的。
(2)对于固定的和b,以及相同的迭代误差条件,两种迭代方法的迭代次数不同。高斯—赛德尔迭代法对应的迭代次数要比雅可比迭代法对应的迭代次数少。这是因为高斯—赛德尔迭代法对应的收敛更快。
(3)对于同一个b,不同的对应的迭代次数不同,且距离解x越远,所需的迭代次数越多。但最后得到的解是一样的。而对于不同的b,最后得到的解一般是不同的。
(3)对于同一个,相同b,给定不同的迭代误差要求,最后得到的解可能略有偏差。误差精度越高,所需的迭代次数越多。误差条件变严格后需要更多的迭代次数以使解满足当前的误差条件。
(5)通过数据表和图像还可以观察到,两种迭代法的收敛情况类似。迭代前几步变化较大,后面均在小范围内变化,与真实解已十分接近,直到满足相应的迭代误差条件时停止迭代过程。
(6)随着矩阵A对角元素的增大,在b和相同以及迭代误差条件相同的前提下,迭代次数越来越少,即收敛速度越来越快。这是由与之对应的的值越来越小所决定的。
值更小,
故与之对应的序列
【题目2】(课本习题第五章第5题)
设国民经济由农业、制造业和服务业三个部门构成,已知某年它们之间的投入产
出关系、外部需求、初始投入等如表7所示。
表7
(1) 如果今年对农业、制造业和服务业的外部需求分别为50、150、100亿元,
问这部门的总产出分别应为多少?
(2) 如果三个部门的外部需求分别增加一个单位,问它们的总产出分别增加多
少?
【第(1)问求解】
由表7得到三个部门的投入产出表如表8所示:
表8
xi
设有n个部门,记一定时期内第i个部门的总产出为投入为
xij
,其中对第j个部门的
,外部需求为
di
,则
xixijdi,i1,2,...,n.
j1n
根据投入系数
aij
的定义,有
xijaijxj,i,j1,2,...,n,
即
aij
是第j个部门的单位产出所需的第i个部门的投入。由此可得
xiaijxjdi,i1,2,...,n.
j1n
,需求向量
记投入系数矩阵
A(aij)nn
,产出向量
x(x1,x2,...,xn)T
d(d1,d2,...,dn)T
,则上式可写作xAxd,或
(IA)xd.
已知今年对农业、制造业和服务业的外部需求分别为50亿元,150亿元,100亿元,求这三个部门的总产出: 【matlab程序】
%------------------------------作业题5_5脚本M文件源程序---------------------------- clear all;clc;
A=[0.15 0.1 0.2;0.3 0.05 0.3;0.2 0.3 0]; % 设定系数矩阵A d=[50 150 100]'; % 给定外部需求 b=eye(3)-A;
x=b\d % 求解方程
dx=inv(b)
得到结果如下:
x 139.2801, 267.6056, 208.1377
T
1.34590.25040.3443
dx0.56341.26760.4930
0.43820.43041.2167
即三个部门的总产出应分别为:139.2801、267.6056、208.1377(亿元)。
【第(2)问求解】
当三个部门的外部需求分别增加1个单位时,求它们总产出相应增加的量: 由(IA)xd可得
x(IA)1d.
表明总产出x对外部需求d是线性的,所以当d增加1个单位(记作d)时,x
1
x(IA)d. 的增量为
于是从dx=inv(b)矩阵中便可得三个部门的外部需求分别增加一个单位时,各部
门总产出分别增加量。
1.34590.25040.3443
dx0.56341.26760.4930
0.43820.43041.2167
结果如下:
T1
d(1,0,0)(IA)x若农业的外部需求增加1个单位,,为的第1列
x=1.3459,0.5634,0.4382;即当农业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加1.3459, 0.5634, 0.4382单位。
若制造业的外部需求增加1个单位,d(0,1,0),x为(IA)的第2列
T
1
T
x=0.2504,1.2676,0.4304;即当制造业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加0.2504,1.2676,0.4304单位。
若服务业的外部需求增加1个单位,d(0,0,1),x为(IA)的第3列
T
1
T
x=0.3443,0.4930,1.2167;即当服务业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加0.3443,0.4930,1.2167单位。
T
【本题小结】
本题是对实际问题的数学建模,模型较简单,问题最后转化为线性方程组求解问题。因为方程并不复杂,我采用直接法求解。很快得到了答案。
【题目3】(课本习题第五章第7题)
输电网络:一种大型输电网络可简化为图5.5所示电路,其中R1,R2,…,Rn表示负载电阻,r1,r2,…,rn表示线路内阻,I1,I2,…,In表示负载上的电流。设电源电压为V。
(1)列出求各负载电流I1,I2,…,In的方程;
(2)设R1=R2=…=Rn=R,r1=r2=…=rn=r,在r=1,R=6,V=18,V=18,n=10的情况下求I1,I2,…,In及总电流I0。
图5.5
【第(1)问求解】
问题分析及模型建立:
IIiIi1Inr
通过对以上电路的分析可得,流过电阻i的电流ri。对于每一个
由Ri,ri+1,Ri1组成的闭合回路,根据基尔霍夫回路电压法,有
-IiRi+Iri+1ri+1Ii1Ri10。对于由电源V,R1,r1组成的回路,有
-V+I1R1Ir1r10
。
因此,可列方程如下:
I1R1(I1+I2+In)r1V0
-IRIR(III)r0112223n2
-I2R2I3R3(I3I4In)r30 -In-1Rn1InRnInrn0
整理后可得,
R11(1+)I+I+IV12nr1r1
R2R1
-I(1)I2I3In=0r1
r22
RR-2I(13)III=0
34n
r32r3
-Rn1I(1Rn)I0
n-1n
rn
rn变成矩阵的形式,即AXb,其中,
1+R1 1 1 … … 1 r1
R1R2r 12r 1 … 12 A= R
2r 3
… Rn1r 1RnnrnIV1
I
2rX
1I03
b0
,。 In-10In
最后列出求各负载电流I1,I2,…,In的方程如下:
R1I1+ 1 1 … … 1 r1V
1
r1
R2R1
r 1r 1 … 1 I20
22
R2
I30 r3
… In-10
Rn1Rn
1 In0 rnrn
【第(2)问求解】
由题意可得, R1=R2=…=Rn=R,r1=r2=…=rn=r,r=1,R=6,V=18, n=10因此,
矩阵A变作
7 1 1 … 1
06 7 1 … 1
6 矩阵b变作 b=0 A=
0
6 7
7 1 1 … 1 I1186 7 1 … 1 I02
整理方程得到 6 I3=0
In-10 6 7In0
18
其中n=10。
【Matlab程序】
%------------------------------作业题5_7脚本M文件源程序---------------------------- clear all;clc; n=10;
A1=sparse(1:n,1:n,7,n,n); A2=sparse(2:n,1:n-1,-6,n,n); A3=triu(ones(10,10),1);
A=A1+A2+A3; % 设定系数矩阵A b=[18,zeros(1,9)]'; % 设定列向量b x=[A\b]' % 求解各支路电流值 I0=sum(x) % 求解电流I0
得到如下结果:
表9:r=1,R=6,V=18,V=18,n=10时I1,I2,…,In,I0的值
所以要求的I1,I2,…,In及总电流I0值如表9所示。
【本题小结】
本题是电路问题的求解。可以采用差分方程求解,也可以用线性方程组求解。我采用线性方程组求解得到答案。因为问题中矩阵A并不复杂,因此采用直接法求解方程组,模型建好后很快便得到了答案。
【实验心得、体会】
本次实验作为数学实验课程的第三次作业,让我们掌握了线性方程组的数值解法相关基础知识,同时我进一步了解了matlab。本次实验第一个题目难度较大其余两题目都比较轻松。通过第一题,我深刻理解了雅克比和高斯——赛德尔迭代法,增长了许多知识。实验中最令我困惑的是为什么很简单的一个矩阵求解(线性方程组求解)问题要用雅克比和高斯——赛德尔迭代法求解,比如第一问中已知矩阵A与b,直接用x=A\b便可得解,我们却要绕一个大弯得到x的解,后来我逐渐明白了,直接法适用于求解中小型(n
我尝试着对每个问题用学到的不同的方法去解决(如对线性方程组分别用直接法中的高斯消元法和LU分解以进行求解,并和两种迭代方法进行比较),同时尝试着从不同的方面和角度去分析,并尽可能地发掘出更多值得去思考和分析的东西,在这个过程中不仅巩固了此次所学到的了线性方程组的数值解等相关知识,也不仅进一步熟练了Matlab的操作,更重要的是通过这次作业进行了一次思维上的极大锻炼,进一步熟悉了通过数学建模思想和数学软件解决一些和现实密切相关的问题。
总之,本次实验收获挺大的。