线性方程组的数值解法

《大学数学实验》作业

线性方程组的数值解法

班级: 姓名: 学号: 日期:

目录

目录 .................................................................................................................................................. 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,其中AR

2020

,定义为

1/21/43

1/231/21/4

1/41/231/21/4

1/41/23

A1/4











31/21/4

1/231/21/41/23

试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数

矩阵性质对收敛速度的影响。实验要求:

(1)选取不同初始向量x和不同的方程组右端向量b,给定迭代误差要求,用雅(0)

克比和高斯——赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;

(2)取定右端向量b和初始向量x(0),将A的主对角元素成倍增长若干次,非主

(k1)(k)5|xx||10对角元素不变,每次用雅克比迭代法,要求迭代误差满足,

比较收敛速度,分析现象并得出你的结论。

【第(1)问求解】

【模型分析】

对于矩阵A,将其分解为:

A=D-L-U,

其中:

30D=

00

0

000

30000

000, L=

030

0

0030

0

000

,U=

00



00

000

000

0

0



0

00 0

其中D、L、U均为20×20的方阵。

则雅可比迭代公式为

x(k1)D1(LU)x(k)D1b, k0,1,2...

高斯—赛德尔迭代公式为:

x(k1)D1(Lx(k1)Ux(k))D1b, k0,1,2...

选取初始向量x

(0)

=[0,0,00],b=[1,1,11],迭代误差为:x(k1)x(k)

105,

用雅可比迭代公法和高斯—赛德尔迭代法对方程组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,00],b=[1,1,11]雅可比迭代公式计算部分结果

用高斯—赛德尔迭代法计算,共迭代11次,得到的部分迭代结果如下表2:

表2:取x

(0)

=[0,0,00],b=[1,1,11]高斯—赛德尔迭代法计算部分结果

x(0)=[0,0,00],b=[1,1,11],x(k1)x(k)

105时,为了更清晰比较两种迭代

方法收敛快慢,作出解x1随迭代次数变化的图像:

x1在两种迭代过程中的变化图像

x1

迭代次数

图1: 解x1随迭代次数变化图像

【结果分析】

由以上数据图表可知,x

(0)

=[0,0,00],b=[1,1,11],x(k1)x(k)

105时,解

向量序列均是收敛的,都收敛到精确解,不过收敛速度不同。雅克比迭代法共迭代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=D1(LU)。在Matlab中运行eig(D_inv*(L+U)语句,得到B=D1(LU)的全部特征值为-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=D1(LU)。利用Matlab的norm命令可以求得B的∞-范数为0.5000。对于高斯—赛德尔迭代法,

。同样可以求得B的

∞-范数为0.2802。即高斯—赛德尔迭代法对应的序列

值更小,故与之对应的

收敛更快,从而对于一定的迭代误差条件,高斯—赛德尔迭代法比雅可

比迭代法法的迭代次数要少。

可见,上述分析与之前的运行结果是相一致的。

【拓展对比】

本题采用的是迭代法求解线性方程组。也可以通过直接法求解线性方程组,

如高斯消元法,LU分解等。可以在Matlab中利用x=A\b或者利用递推关系式进行求解。求出来的结果与用迭代法求出来的解略有出入。这是由求解过程中的误差所致。

【第(2)问求解】

取定方程组右端向量b=[1,2,320],同时取定初始向量x(0)=[0,0,00],并将迭代误差设为为:x(k1)x(k)

105,将A的主对角线元素分别增长为6,

12,24,48,非主对角元素不变,每次用雅可比迭代法进行迭代计算,观察收敛速度的变化,得到结果如下面的表6:

表6: A主对角元素加倍时收敛速度的变化(b=

,x(0)=

【结果分析】

由表6可以明显看到,随着矩阵A对角元素的增大,对于雅可比迭代法,与之对应的

的值越来越小。在b和

相同以及迭代误差条件相同的前提

下,迭代次数越来越少,即收敛速度越来越快。

分析: 由前一问对收敛速度的分析可得, 由于矩阵的谱半径不超过它的任一种范数,所以若收敛,且有误差估计式

,则迭代公式

可见,

越小,序列

收敛越快。

对于雅可比迭代法,B=D1(LU)。利用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

,则

xixijdi,i1,2,...,n.

j1n

根据投入系数

aij

的定义,有

xijaijxj,i,j1,2,...,n,

aij

是第j个部门的单位产出所需的第i个部门的投入。由此可得

xiaijxjdi,i1,2,...,n.

j1n

,需求向量

记投入系数矩阵

A(aij)nn

,产出向量

x(x1,x2,...,xn)T

d(d1,d2,...,dn)T

,则上式可写作xAxd,或

(IA)xd.

已知今年对农业、制造业和服务业的外部需求分别为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

dx0.56341.26760.4930

0.43820.43041.2167

即三个部门的总产出应分别为:139.2801、267.6056、208.1377(亿元)。

【第(2)问求解】

当三个部门的外部需求分别增加1个单位时,求它们总产出相应增加的量: 由(IA)xd可得

x(IA)1d.

表明总产出x对外部需求d是线性的,所以当d增加1个单位(记作d)时,x

1

x(IA)d. 的增量为

于是从dx=inv(b)矩阵中便可得三个部门的外部需求分别增加一个单位时,各部

门总产出分别增加量。

1.34590.25040.3443

dx0.56341.26760.4930

0.43820.43041.2167

结果如下:

T1

d(1,0,0)(IA)x若农业的外部需求增加1个单位,,为的第1列

x=1.3459,0.5634,0.4382;即当农业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加1.3459, 0.5634, 0.4382单位。

若制造业的外部需求增加1个单位,d(0,1,0),x为(IA)的第2列

T

1

T

x=0.2504,1.2676,0.4304;即当制造业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加0.2504,1.2676,0.4304单位。

若服务业的外部需求增加1个单位,d(0,0,1),x为(IA)的第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)问求解】

问题分析及模型建立:

IIiIi1Inr

通过对以上电路的分析可得,流过电阻i的电流ri。对于每一个

由Ri,ri+1,Ri1组成的闭合回路,根据基尔霍夫回路电压法,有

-IiRi+Iri+1ri+1Ii1Ri10。对于由电源V,R1,r1组成的回路,有

-V+I1R1Ir1r10

因此,可列方程如下:

I1R1(I1+I2+In)r1V0

-IRIR(III)r0112223n2

-I2R2I3R3(I3I4In)r30 -In-1Rn1InRnInrn0

整理后可得,

R11(1+)I+I+IV12nr1r1

R2R1

-I(1)I2I3In=0r1

r22

RR-2I(13)III=0

34n

r32r3

-Rn1I(1Rn)I0

n-1n

rn

rn变成矩阵的形式,即AXb,其中,

1+R1 1 1 … … 1 r1

R1R2r 12r 1 …  12 A= R

2r  3

  …     Rn1r 1RnnrnIV1

I

2rX

1I03

b0

,。 In-10In

最后列出求各负载电流I1,I2,…,In的方程如下:







R1I1+ 1 1 … … 1 r1V

1

r1

R2R1

r 1r 1 …  1 I20

22

R2

  I30 r3

  …   In-10 

Rn1Rn

1 In0 rnrn

【第(2)问求解】

由题意可得, R1=R2=…=Rn=R,r1=r2=…=rn=r,r=1,R=6,V=18, n=10因此,

矩阵A变作

7 1 1 … 1 

06 7 1 … 1 



 6   矩阵b变作 b=0 A=    

0

 6 7

7 1 1 … 1 I1186 7 1 … 1 I02

整理方程得到  6  I3=0

    In-10 6 7In0

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,其中AR

2020

,定义为

1/21/43

1/231/21/4

1/41/231/21/4

1/41/23

A1/4











31/21/4

1/231/21/41/23

试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数

矩阵性质对收敛速度的影响。实验要求:

(1)选取不同初始向量x和不同的方程组右端向量b,给定迭代误差要求,用雅(0)

克比和高斯——赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;

(2)取定右端向量b和初始向量x(0),将A的主对角元素成倍增长若干次,非主

(k1)(k)5|xx||10对角元素不变,每次用雅克比迭代法,要求迭代误差满足,

比较收敛速度,分析现象并得出你的结论。

【第(1)问求解】

【模型分析】

对于矩阵A,将其分解为:

A=D-L-U,

其中:

30D=

00

0

000

30000

000, L=

030

0

0030

0

000

,U=

00



00

000

000

0

0



0

00 0

其中D、L、U均为20×20的方阵。

则雅可比迭代公式为

x(k1)D1(LU)x(k)D1b, k0,1,2...

高斯—赛德尔迭代公式为:

x(k1)D1(Lx(k1)Ux(k))D1b, k0,1,2...

选取初始向量x

(0)

=[0,0,00],b=[1,1,11],迭代误差为:x(k1)x(k)

105,

用雅可比迭代公法和高斯—赛德尔迭代法对方程组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,00],b=[1,1,11]雅可比迭代公式计算部分结果

用高斯—赛德尔迭代法计算,共迭代11次,得到的部分迭代结果如下表2:

表2:取x

(0)

=[0,0,00],b=[1,1,11]高斯—赛德尔迭代法计算部分结果

x(0)=[0,0,00],b=[1,1,11],x(k1)x(k)

105时,为了更清晰比较两种迭代

方法收敛快慢,作出解x1随迭代次数变化的图像:

x1在两种迭代过程中的变化图像

x1

迭代次数

图1: 解x1随迭代次数变化图像

【结果分析】

由以上数据图表可知,x

(0)

=[0,0,00],b=[1,1,11],x(k1)x(k)

105时,解

向量序列均是收敛的,都收敛到精确解,不过收敛速度不同。雅克比迭代法共迭代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=D1(LU)。在Matlab中运行eig(D_inv*(L+U)语句,得到B=D1(LU)的全部特征值为-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=D1(LU)。利用Matlab的norm命令可以求得B的∞-范数为0.5000。对于高斯—赛德尔迭代法,

。同样可以求得B的

∞-范数为0.2802。即高斯—赛德尔迭代法对应的序列

值更小,故与之对应的

收敛更快,从而对于一定的迭代误差条件,高斯—赛德尔迭代法比雅可

比迭代法法的迭代次数要少。

可见,上述分析与之前的运行结果是相一致的。

【拓展对比】

本题采用的是迭代法求解线性方程组。也可以通过直接法求解线性方程组,

如高斯消元法,LU分解等。可以在Matlab中利用x=A\b或者利用递推关系式进行求解。求出来的结果与用迭代法求出来的解略有出入。这是由求解过程中的误差所致。

【第(2)问求解】

取定方程组右端向量b=[1,2,320],同时取定初始向量x(0)=[0,0,00],并将迭代误差设为为:x(k1)x(k)

105,将A的主对角线元素分别增长为6,

12,24,48,非主对角元素不变,每次用雅可比迭代法进行迭代计算,观察收敛速度的变化,得到结果如下面的表6:

表6: A主对角元素加倍时收敛速度的变化(b=

,x(0)=

【结果分析】

由表6可以明显看到,随着矩阵A对角元素的增大,对于雅可比迭代法,与之对应的

的值越来越小。在b和

相同以及迭代误差条件相同的前提

下,迭代次数越来越少,即收敛速度越来越快。

分析: 由前一问对收敛速度的分析可得, 由于矩阵的谱半径不超过它的任一种范数,所以若收敛,且有误差估计式

,则迭代公式

可见,

越小,序列

收敛越快。

对于雅可比迭代法,B=D1(LU)。利用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

,则

xixijdi,i1,2,...,n.

j1n

根据投入系数

aij

的定义,有

xijaijxj,i,j1,2,...,n,

aij

是第j个部门的单位产出所需的第i个部门的投入。由此可得

xiaijxjdi,i1,2,...,n.

j1n

,需求向量

记投入系数矩阵

A(aij)nn

,产出向量

x(x1,x2,...,xn)T

d(d1,d2,...,dn)T

,则上式可写作xAxd,或

(IA)xd.

已知今年对农业、制造业和服务业的外部需求分别为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

dx0.56341.26760.4930

0.43820.43041.2167

即三个部门的总产出应分别为:139.2801、267.6056、208.1377(亿元)。

【第(2)问求解】

当三个部门的外部需求分别增加1个单位时,求它们总产出相应增加的量: 由(IA)xd可得

x(IA)1d.

表明总产出x对外部需求d是线性的,所以当d增加1个单位(记作d)时,x

1

x(IA)d. 的增量为

于是从dx=inv(b)矩阵中便可得三个部门的外部需求分别增加一个单位时,各部

门总产出分别增加量。

1.34590.25040.3443

dx0.56341.26760.4930

0.43820.43041.2167

结果如下:

T1

d(1,0,0)(IA)x若农业的外部需求增加1个单位,,为的第1列

x=1.3459,0.5634,0.4382;即当农业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加1.3459, 0.5634, 0.4382单位。

若制造业的外部需求增加1个单位,d(0,1,0),x为(IA)的第2列

T

1

T

x=0.2504,1.2676,0.4304;即当制造业的需求增加1个单位时,农业、制造业和服务业的总产出应分别增加0.2504,1.2676,0.4304单位。

若服务业的外部需求增加1个单位,d(0,0,1),x为(IA)的第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)问求解】

问题分析及模型建立:

IIiIi1Inr

通过对以上电路的分析可得,流过电阻i的电流ri。对于每一个

由Ri,ri+1,Ri1组成的闭合回路,根据基尔霍夫回路电压法,有

-IiRi+Iri+1ri+1Ii1Ri10。对于由电源V,R1,r1组成的回路,有

-V+I1R1Ir1r10

因此,可列方程如下:

I1R1(I1+I2+In)r1V0

-IRIR(III)r0112223n2

-I2R2I3R3(I3I4In)r30 -In-1Rn1InRnInrn0

整理后可得,

R11(1+)I+I+IV12nr1r1

R2R1

-I(1)I2I3In=0r1

r22

RR-2I(13)III=0

34n

r32r3

-Rn1I(1Rn)I0

n-1n

rn

rn变成矩阵的形式,即AXb,其中,

1+R1 1 1 … … 1 r1

R1R2r 12r 1 …  12 A= R

2r  3

  …     Rn1r 1RnnrnIV1

I

2rX

1I03

b0

,。 In-10In

最后列出求各负载电流I1,I2,…,In的方程如下:







R1I1+ 1 1 … … 1 r1V

1

r1

R2R1

r 1r 1 …  1 I20

22

R2

  I30 r3

  …   In-10 

Rn1Rn

1 In0 rnrn

【第(2)问求解】

由题意可得, R1=R2=…=Rn=R,r1=r2=…=rn=r,r=1,R=6,V=18, n=10因此,

矩阵A变作

7 1 1 … 1 

06 7 1 … 1 



 6   矩阵b变作 b=0 A=    

0

 6 7

7 1 1 … 1 I1186 7 1 … 1 I02

整理方程得到  6  I3=0

    In-10 6 7In0

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的操作,更重要的是通过这次作业进行了一次思维上的极大锻炼,进一步熟悉了通过数学建模思想和数学软件解决一些和现实密切相关的问题。

总之,本次实验收获挺大的。


相关文章

  • 非线性方程和线性方程的数值解
  • 非线性方程的几种数值解法 郭摇 摘要:现在许多科学理论和工程技术问题都最终要化成非线性方程f(x)=0或非线性方程组F(x)=0的求解.对于非线性方程的数值解法有很多种,如:二分法,简易牛顿法,斯蒂芬森迭代法,牛顿迭代法,割线法等. 关键字 ...查看


  • 数值分析学习心得体会
  • 数值分析学习感想 一个学期的数值分析,在老师的带领下,让我对这门课程有了深刻的理解和感悟.这门 课程是一个十分重视算法和原理的学科,同时它能够将人的思维引入数学思考的模式,在处 理问题的时候,可以合理适当的提出方案和假设.他的内容贴近实际, ...查看


  • 常微分方程常用数值解法
  • 第一章 绪论 1.1 引言 常微分方程是现代数学的一个重要分支,是人们解决各种实际问题的有效工具.微分方程的理论和方法从17世纪末开始发展起来,很快成了研究自然现象的强有力工具,在17到18世纪,在力学.天文.科学技术.物理中,就已借助微分 ...查看


  • 线性方程组的数值解法 1
  • [实习目的] 1 通过实习进一步掌握高斯消去法,列主元高斯消去法的基本思想: 2 通过实习进一步掌握高斯消去法,列主元高斯消去法的计算步骤,并能灵活应用: 3 通过对高斯消去法,列主元高斯消去法的调试练习,进一步体会各种算法的特点: 4 通 ...查看


  • 高斯列主元消去法解线性方程组
  • 线性方程组的数值解法 --高斯列主元消去法解线性方程组的MATLAB实现 班级:MATH 20XX 学号:xxxxxxxxxx 姓名:LI Wen 高斯列主元消去法解线性方程组的MATLAB实现 摘 要 在自然科学和工程技术中许多问题的解决 ...查看


  • 数值分析在大地测量学与测量工程中的应用
  • 数值分析在大地测量中的应用 摘要:数值分析是一门研究最基本最通用的数值算法,即研究基本数学问题的适合计算机求解的数值方法.数值分析,在测绘领域取得了重要而广泛的应用.本文简要阐述了数值分析的基本特点.发展规律和发展趋势,结合实例,主要对数值 ...查看


  • 常微分方程差分解法.入门.多解法
  • 毕业论文 题 目 抛物型方程的差分解法 学 院 专 业 信息与计算科学 班 级 学 生 王丹丹 学 号 指导教师 王宣欣 二〇一二 年 五 月二十五日 摘 要 偏微分方程的数值解法在数值分析中占有重要的地位,很多科学技术问题的数值计算包括了 ...查看


  • 线性方程组的解法
  • 线性方程组的解法 1 引言 在科学研究和大型工程设计中出现了越来越多的数学问题,而这些问题往往需要求数值解.在进行数值求解时,经离散后,常常归结为求解形如Ax= b的大型线性方程组.而如插值公式,拟合公式等的建立,微分方程差分格式的构造等, ...查看


  • 算法大全第15章常微分方程的解法
  • 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步,通常需要求出方程的解来说明实际现象,并加以检验.如果能得到解析形式的解固然是便于分析和应用的,但是我们知道,只有线性常系数微分方程,并且自由项是某些特殊类型的函数时,才可以肯定得 ...查看


  • 斜拉索非线性振动的奇异摄动解法
  • 西 南 交 通 大 学 学 报第41卷 第3期Vol . 41 No . 3 2006 年6月Jun . 2006JOURNAL OF S OUT HW EST J I A OT ONG UN I V ERSI TY 文章编号:025822 ...查看


热门内容