max=abs（A（j,i））;
m=j;
end
end
if（m~=i）
fork=i:
n
c（k）=A（i,k）;
A（i,k）=A（m,k）;
A（m,k）=c（k）;
end
d1=B（i）;
B（i）=B（m）;
B（m）=d1;
end
fork=i+1:
n
forj=i+1:
n
A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;
end
B（k）=B（k）-B（i）*A（k,i）/A（i,i）;
A（k,i）=0;
end
end
X（n）=B（n）/A（n,n）;（忘记点标点“；”）
fori=n-1:
-1:
1
sum=0;
forj=i+1:
n
sum=sum+A（i,j）*X（j）;
end
X（i）=（B（i）-sum）/A（i,i）;
end
X
计算结果为：
X=（0.999999；1.000000；1.000002；1.000000）
【命令】：
inv（A）
ans=
-1100
1-210
01-21
001-2
第二题：
选用一种编程语言实现对一般三对角线性方程组的追赶法。
具体算例可以参考书P61.第2题，也可以自己设计算例。
1.对追赶法的叙述：
对于系数矩阵是稀疏的三对角矩阵的方程组，追赶法是用来求解之的专用方法。
对于三对角方程组，追赶法比Gauss消去法的计算量要小的多。
如果A满足三角分解的条件，可将A分解为LU，分解后将线性方程组化为求解方程组Ly=d和Ux=y，从而得到三对角方程组的解。
追赶法是由于矩阵中出现了大量零元素，从而使得计算公式简化，计算量大为减少，其乘除法的运算量为5n-4.
2.追赶法的算法：
（1）输入a=（a2,...,an）,b=（b1,b2,...,bn）,c=（c1,c2,...,cn-1）,d=（d1,d2,...,dn）,n.
（2）对于i=2,3,...，n,置ai/bi-1=li→ai,bi-ci-1ai=ui→bi,di-aidi-1=yi→di.
（3）置dn/bn=xn→dn
（4）对于i=n-1,n-2,...,1,置（di-cidi+1）/bi=xi→di
（5）输入d=x,停机.
注：
因为三对角矩阵的非零元素都集中在三条对角线上，因此只用三个向量a,b,c来存储系数矩阵，这样可以把二维数据的存储变为一维数据存储。
此外，因为追赶法的LU分解是特殊的分解，它只有二个向量数据p,q。
与Doolittle分解一样考虑，为节省计算机存储单元，用原来系数矩阵A的向量a,b来贮存p,q，不但省存储单元，而且可以减少计算量。
3.具体算例：
试用MATLAB软件编程实现追赶法求解三对角方程组的算法，并考虑如下梯形电阻电路问题：
其中电路中的各个电流{i1,i2,…，i8}满足下列线性方程组：
2i1-2i2 =V/R
-2i1+5i2-2i3=0
-2i2+5i3-2i4=0
-2i3+5i4-2i5=0
-2i4+5i5-2i6=0
-2i5+5i6-2i7=0
-2i6+5i7-2i8=0
-2i7+5i8=0
设V=220V，R=27Ω，求各段电路的电流量。
数值结果i1=149.0717i2=70.4618i3=31.1568i4=12.8624i5=4.9497i6=1.7169i7=0.4772i8=0.0477
4.总结
由于追赶法来源于Gauss法，也存在其缺点即对角元素为0时无法计算。
在用matlab编程时经常出错，所以一定要有耐心。
但追赶法由于矩阵中出现了大量零元素，计算简单，而且其乘除法的运算量小.特殊的求解过程，可以节省存储空间。
5.程序
Clear
m=8;
n=8;
V=220;
R=27;
A=zeros（m,n）;
c=zeros（n,1）;
i=zeros（n,1）;
A（1,1）=2;
fork=1:
n-1
A（k+1,k）=-2;
A（k+1,k+1）=5;
A（k,k+1）=-2;
end
c
（1）=V/R;
fork=1:
n-1
A（k+1,k）=A（k+1,k）/A（k,k）;
A（k+1,k+1）=A（k+1,k+1）-A（k+1,k）*A（k,k+1）;
end
fork=2:
n
c（k）=c（k）-A（k,k-1）*c（k-1）;
end
i（n）=c（n）/A（n,n）;
fork=n-1:
-1:
1
i（k）=（c（k）-A（k,k+1）*i（k+1））/A（k,k）;
end
fork=1:
8
fprintf（'i（%d）=',k）;
fprintf（'%9.8f',i（k））;
end
程序运行结果：
i1=8.14777517i2=4.07370109i3=2.03647757i4=1.01749282i5=0.50725449i6=0.25064339i7=0.11935400i8=0.04774160
第三题：
选用一种编程语言实现对一般线性方程组得Jaocobi迭代、Gauss-Seidel迭代和SOR迭代算法，设计算例从数值上比较三种迭代法的收敛结果与收敛速度，并从理论上对所得结论加以分析。
算例可以参考P89.1.
1.对算法的叙述
当解答高阶方程组时，如果用直接方法求解，即使系数矩阵是稀疏的，其稀疏性也很难保持。
迭代法由于能保持矩阵的稀疏性，计算简单，编程容易等优点，并在有些情况下收敛较快，从而成为求解高阶方程组的方法。
Jaocobi迭代法中要求系数矩阵A非奇异，且aii≠0。
迭代公式为x（k+1）=Bx（k）+g，迭代产生向量序列｛x（k）｝。
在Jaocobi迭代法中，如果用新的分量xl（k+1）代替旧的分量xl（k），这就是Gauss-Seidel迭代法。
为了加速迭代向量序列｛x（k）｝的收敛速度，在Gauss-Seidel迭代法的基础上，引入参数ω，可以得到松弛法。
2.对算法的具体描述
（1）Jacobi迭代法具体步骤：
线性方程组的系数矩阵A如果可逆且主对角元素
均不为零,将对角元素组成矩阵
将方程变形，且将x的系数化为1，即变为下面的形式：
其中
.
以
为迭代矩阵的迭代法（公式）
称为雅可比（Jacobi）迭代法（公式）,用向量的分量来表示。
（2）Gauss-Seidel迭代法具体步骤：
把矩阵A分解成
其中
分别为
的主对角元除外的下三角和上三角部分,于是,方程组便可以写成
其中
以
为迭代矩阵构成的迭代法（公式）
称为高斯—赛德尔迭代法（公式）。
（3）SOR迭代算法的具体步骤：
把Gauss-Seidel迭代
的结果作为中间值，记为
。
SOR方法是将
与上次计算的结果
做加权平均作为最后结果。
迭代格式为：
算法：
1．
2．当
时，
，结果仍然存储在
中。
迭代次数
3．计算误差
（真解已知）
4．如果
，则已达到精确度要求，否则继续第2步
3.具体算例：
A=[1,9,-1,2;2,-1,7,3;10,1,2,3;3,2,3,12]b=[-27;14;12;-17]初始向量x=[0;0;0;0];
Jacobi迭代法结果X=[1.0000；-2.0000；3.0000；-2.0000]
Gauss-Seidel迭代法结果X=[1.0000；-2.0000；3.0000；-2.0000]
共轭梯度法结果X=[1.0000；-2.0000；3.0000；-2.0000]
4.总结
从以上结果可以看出在求解相同问题时，sor法的迭代次数要比Jacobi迭代法少很多，计算量小许多。
此外可以看出松弛因子w的选取对迭代次数的影响也十分大。
在实际计算时，最优松弛因子很难事先确定，一般可用试算法取近似最优值。
高斯—塞德尔迭代法比雅可比迭代法收敛快（达到同样的精度所需迭代次数少）,但这个结论,在一定条件下才是对的,甚至有这样的方程组,雅可比方法收敛，而高斯—塞德尔迭代法却是发散的。
在计算机工具的帮助下，可以很快的实现方程的求解。
但在编程的时候可能一个标点的疏忽就可能导致错误，一定要细心。
5.程序
设计程序如下：
（1）Jacobi：
b=[-27;14;12;-17]
x=[0;0;0;0]
k=0;
r=1;
e=0.000001
A=[1,9,-1,2;2,-1,7,3;10,1,2,3;3,2,3,12]
D=diag（diag（A））;（矩阵的维数弄错）
B=inv（D）*（D-A）;
f=inv（D）*b;
p=max（abs（eig（B）））;
ifp>=1
'迭代法不收敛'
return
end
whiler>e
x0=x;
x=B*x0+f;
k=k+1;
r=norm（x-x0,inf）;
end
x
k
计算结果：
x=
1.0000
-2.0000
3.0000
-2.0000
k=65
（2）Gauss-Seidel迭代
设计程序如下：
A=[1,9,-1,2;2,-1,7,3;10,1,2,3;3,2,3,12]
x=[0;0;0;0;0];
b=[-27;14;12;-17]
c=0.000001
L=-tril（A,-1）
U=-triu（A,1）
D=（diag（diag（A）））
X=inv（D-L）*U*x+inv（D-L）*b;
k=1;
whilenorm（X-x,inf）>=c
x=X;
X=inv（D-L）*U*x+inv（D-L）*b;
k=k+1;
end
X
k
计算结果：
X=
1.0000
-2.0000
3.0000
-2.0000
k=37
（3）SOR法：
A=[1,9,-1,2;2,-1,7,3;10,1,2,3;3,2,3,12]x=[0;0;0;0;0]b=[-27;14;12;-17];
e=0.000001
w=1.44;
L=-tril（A,-1）
U=-triu（A,1）
D=（diag（diag（A）））
X=inv（D-w*L）*（（1-w）*D+w*U）*x+w*inv（D-w*L）*b
n=1;
whilenorm（X-x,inf）>=e
x=X;
X=inv（D-w*L）*（（1-w）*D+w*U）*x+w*inv（D-w*L）*b;
n=n+1;
end
X
n
计算结果：
X=
1.0000
-2.0000
3.0000
-2.0000
n=22
结果分析：
由迭代次数的比较可知该情况下Jacobi迭代法比Gauss-Seidel迭代法收敛的慢，SOR迭代法收敛的快，但松弛因子难确定。
第四题、选用一种编程语言实现幂法和反幂法。
自行设计算例，也可参看P122.4.
1.对幂法和反幂法的叙述
幂法主要用于求矩阵的按模最大的特征值与相应的特征向量，它是通过迭代产生向量序列，由此计算特征值和特征向量。
幂法的收敛速度是线性的，而且依赖于收敛因子，收敛因子接近1时很慢，可以进行加速。
反幂法是计算矩阵按模最小的特征值及特征向量的方法，也是修正特征值、求相应特征向量的最有效的方法。
2.具体算法过程
幂法的算法：
（1）输入A=（aij）,初始向量x=（x1,x2,...,xn）,误差限ε，最大迭代次数N；
（2）置k=1，μ=0；
（3）求整数r，使得|xr|=max,xr→α;
（4）计算y=x/α,x=Ay,置xr→λ;
（5）若|λ-μ|<ε，输出λ，x,停机；否则，转6；
（6）若k反幂法的算法：
（1）输入A=（aij）,近似值λ*，初始向量x=（x1,x2,...,xn）,误差限ε，最大迭代次数N；
（2）置k=1，μ=1；
（3）作三角分解（A-λ*I）=LU；
（4）求整数r，使得|xr|=max,xr→α;
（5）计算y=x/α，Lz=y,Ux=z,置xr→β；
（6）若|1/β-1/μ|<ε，则置λ=λ*+1/β，输入λ，x,停机；否则转7；
（7）若k3.具体算例：
设A=[6-126;-21-324;-12-1251]，取x（0）=（1,1,1）T,先用幂法迭代3次，得到A的按模最大特征值的近似值，取λ*为其整数部分，再用反幂法计算A的按模最大特征值的更精确的近似值，要求误差小于10-10。
数值结果幂法x1=[0.5159；-0.8371；-0.1819]反幂法x2=[-0.8126;-1.2191;-0.4063]
]
4.总结
反幂法用来计算矩阵按模最小的特征值及其特征向量，也可用来计算对应与一个给定近似特征值的特征向量。
它的计算过程需要解线性方程组，一般可以采用直接三角分解法。
幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法,特别是用于大型稀疏矩阵，也是计算矩阵谱半径的有效方法。
但它的收敛速度是线性的，比较慢。
5.程序
（1）幂法
>>A=[6-126;-21-324;-12-1251];
N=100;
ep=1e-10;
>>x=[1;1;1];
>>k=0;
mu=0;
whilek<=N
alpha=max（abs（x））;
y=x/alpha;
x=A*y;
lambda=alpha;
ifabs（lambda-mu）break;
end
k=k+1;mu=lambda;
end
x=x/norm（x）
x=
0.5159
-0.8371
-0.1819
（2）反幂法
反幂法
>>A=[6-126;-21-324;-12-1251];
N=100;
ep=1e-10;
x=[1;1;1];
>>v=x;%按模最小的特征向量v
m=0;
l=0;%按模最小的特征值
for（k=1:
N）
yy=A\v;
m=max（abs（yy））;
v=yy/m;
if（abs（m-l）l=l/m;
return;
end
end
fprintf（'\n矩阵A的按模最小的特征值λ=%f\n\n',l）;
fprintf（'λ对应的特征向量=[%f;%f;%f;%f]\n',v）;
矩阵A的按模最小的特征值λ=-9.002
λ对应的特征向量=[-0.8126;-1.2191;-0.4063]
第五题、选用一种编程语言实现逐次分半梯形算法和Simpson算法。
自行设计算例。
1.对逐次分半梯形算法和simpson算法的叙述
复化求积公式的截断误差随着步长的缩小而减小，由给定的精度可以预先确定步长。
但由于误差估计式中出现高阶导数，使问题变得十分困难。
实际计算时，通常是从某个步长h出发计算近似值，然后在前次分割的基础上，再将每个小区间等分，即将积分区间逐次分半，并以前后两次计算结果之差来估计误差，直到求得满足精度的积分近似值。
这样既缩小了步长，又保留了已有的节点，减少了计算量，这就是逐次分半算法。
2.逐次分半梯形算法和simpson算法的具体描述
（1）i.输入a,b,f（x）,ε；
ii.置m=1,h=（b-a）/2,T0=h[f（a）+f（b）]；
iii.置F=0，对k=1,2,...,2m-1F=F+f（a+（2k-1）h）；
iv.T=T0/2+hF；
V.若|T-T0|<3ε，输出I≈T，停机；否则m+1→m,h/2→h,T→T0,转3.
（2）i.输入a,b,f（x）,ε；
ii.置F1=f（a）+f（b）,F2=f（（a+b）/2）,S0=（b-a）（F1+4F2）/6,m=2,h=（b-a）/4;
iii.置F3=0，对k=1,2,...,2m-1,F3=F3+f（a+（2k-1）h）;
Iv.S=h（F1+2F2+4F3）；
v.若|S-S0|<15ε,则输出I≈S,停机;否则,m+1→m,h/2→h,F2+F3→F2,S→S0,转3.
3.具体算例
dx
数值结果
（1）梯形公式的逐次分半算法T4=1.7272219
（2）Simpson的逐次分半算法S4=1.7183188
4.总结
逐次分半算法将积分区间逐次分半，由前后两次计算之差来估计误差，使其满足精度。
对于求积分而言，用逐次分半算法能更快得到满足精度的解。
在进行积分运算时，要注意表达式应该有意义，因此应该注意积分的上下限是不是对被积函数有意义。
5.程序
（1）梯形公式的逐次分半算法
functionI=tixingstep（fx=ex,a=0,b=1,eps=1/2*10-2）;
h=b-a;
c=（a+b）/2;m=1;h=h/2;
f1=feval（fx,a）+feval（fx,b）;f2=0;
s0=0;s=h*（f1+2f2）;
while（abs（s-s0）>=3*eps）
m=m+1;h=h/2;n=2^m;f2=f2+f2;s0=s;
x=a+h:
h:
b-h;f3=sum（feval（fx,x））;
s=h/2*（f1+2*f2）;
end
I=T=1.7272219;
（2）Simpson公式的逐次分半法
functionI=simpsonstep（fx=ex,a=0,b=1,eps=1/2*10-2）%fx为设定的被积函数
h=b-a;c=（a+b）/2;m=1;h=h/2;
f1=feval（fx,a）+feval（fx,b）;f2=feval（fx,c）;f3=0;%feval（）为求函数fx在各点的值
s0=0;s=h/6*（f1+4*f2）;
While（abs（s-s0）>=15*eps）
m=m+1;h=h/2;n=2^m,f2=f2+f3;s0=s;
x=a+h:
2*h:
b~h;f3=sum（feval（fx,x））;
s=h/3*（f1+2*f2+4*f3）;
end
I=s=1.7183188;
升级会员 