数值分析 - 编程作业

《数值分析》实验报告

第二章：解线性方程组的直接方法

2、试用MATLAB软件编程实现追赶法求解三对角方程组的算法，并考虑梯形电阻电路问题，电路如下：

其中电路中的各个电流{i1，i2，?，i8}须满足下列线性方程组：

2i1?2i2 ?VR

?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?，运用求各段电路的电流量。

解：VR?220?8.1481 27上述方程组可用矩阵表示为：

?2??2??0??0?0??0?0???0?25?2000000?25?2000000?25?2000000?25?2000000?25?2000000?25?20??i1??8.1481??i??0?0?2??????0??i3??0?????0??i4??0? ??0??i5??0?????0??i6??0???2??i7??0?????5?????0??i8???MatLab程序： %赋初值；

a=[0 -2 -2 -2 -2 -2 -2 -2]; b=[2 5 5 5 5 5 5 5];

c=[-2 -2 -2 -2 -2 -2 -2]; d=[8.1481 0 0 0 0 0 0 0]; %三对角方程的追赶法

for i=2:8%“追”的过程； a(i)=a(i)/b(i-1);

b(i)=b(i)-c(i-1)*a(i); d(i)=d(i)-a(i)*d(i-1); end;

d(8)=d(8)/b(8);%“赶”的过程； for i=7:-1:1

d(i)=(d(i)-c(i)*d(i+1))/b(i); end;

x=d; x

程序运行结果：

x =

8.1477 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477

?i1??8.1477??i??4.0737??2????i3??2.0365?????i1.01754? (A)。 即I??????i5??0.5073??????i6??0.2506??i??0.1194??7????0.0477???i8????

1、试分别用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法；（3）共轭梯度法解线性方程组

1234??x1??12??10?1??x???27?9?12?3???2????2?173?5??x3???14? ??????x32312?1?17???4?????4?3?5?115????12???x5???

迭代初始向量取x(0)=(0，0，0，0，0)T。 解：实验步骤及程序、结果

取要求达到的精度??10?8。以下程序中的k均表示迭代次数。

（1）Jacobi迭代法MatLab源程序。

format long

A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];

b=[12,-27,14,-17,12];

x0=[0,0,0,0,0];x1=x0; Nmax=1000; k=0;

for i=1:5

sum=0;

for j=1:5

if j~=i

sum=sum+A(i,j)*x0(j);

end; end;

x1(i)=(b(i)-sum)/A(i,i); end;

while abs(norm(x1-x0,inf))>1e-8 & k

x0=x1; for i=1:5 sum=0;

for j=1:5 if j~=i

sum=sum+A(i,j)*x0(j);

end; end;

x1(i)=(b(i)-sum)/A(i,i); end; k=k+1; end;

x1 k

输出结果： x1 =

1.000000024731627 -2.000000022663167 -1.999999993335823 0.999999970325168

k =88

迭代了88次。

2.999999958446139

（2）Gauss-Seidel迭代法MatLab源程序。

A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15]; b=[12,-27,14,-17,12];

x0=[0,0,0,0,0];x1=x0; Nmax=1000; k=1;sum=0;

for j=2:5

sum=sum+A(1,j)*x0(j); end;

x1(1)=(b(1)-sum)/A(1,1);

for i=2:4

sum=0;

for j=1:(i-1)

sum=sum+A(i,j)*x0(j); end;

x1(i)=b(i)-sum; sum=0;

for j=(i+1):5

sum=sum+A(i,j)*x0(j); end;

x1(i)=(x1(i)-sum)/A(i,i);

end; sum=0;

for j=1:4

sum=sum+A(5,j)*x0(j); end;

x1(5)=(b(5)-sum)/A(5,5);

while abs(norm(x1-x0,inf))>1e-8 & k

x0=x1; sum=0;

for j=2:5

sum=sum+A(1,j)*x0(j); end;

x1(1)=(b(1)-sum)/A(1,1);

for i=2:4

sum=0;

for j=1:(i-1)

sum=sum+A(i,j)*x0(j); end;

x1(i)=b(i)-sum; sum=0;

for j=(i+1):5