temp=a(j+1);
a(j+1)=a(j);
a(j)=temp;
end
end
end
%构造矩阵B
fori=1:
1:
30
A(i,i)=1;
ifi>1
fork=1:
1:
i-1
A(i,k)=a(i-k)
end
end
forj=i+1:
1:
30
A(i,j)=a(j-i);
end
end
%构造向量b
b=a';
%gauss消元解方程组
na=size(A,2);
nb=size(b,1);
C=[A,b];
ifna~=nb
error('A和b必须有一致的维数');
end
%消元过程
fori=1:
1:
na-1
forj=i+1:
1:
na
C(j,:
)=C(j,:
)-C(j,i)/C(i,i)*C(i,:
);
end
end
fori=1:
1:
na
forj=1:
1:
na
ifC(i,j)<10*eps
C(i,j)=0;
end
end
end
A=C(:
1:
na);
b=C(:
na+1);
%回带得方程组的解
b1=b;
fori=na:
-1:
1
x(i)=b1(i)/A(i,i);
ifi>1
forj=i-1:
-1:
1
b1(j)=b1(j)-A(j,i)*x(i);
end
end
end
有X=Gauss(A,a)得方程组的解为:
X=[-1.29053.3290-4.56611.49470.3128-0.65172.3771-1.8050
2.2603-1.38700.12841.37740.88731.6156-0.70630.7630
-1.9829-0.3611-1.28300.52570.0903-2.12871.6322-2.4127
1.5899-1.87370.09613.8652-3.47072.5534]’
(4)、平方根法求解程序如下:
functionx=average(A,b)
n=size(A,1);
l(1,1)=sqrt(A(1,1));
fori=2:
n
l(i,1)=A(i,1)/l(1,1);
end
forj=2:
n
sum1=0;
fork=1:
j-1
sum1=sum1+l(j,k)*l(j,k);
end
l(j,j)=sqrt(A(j,j)-sum1);
fori=j+1:
n
sum2=0;
fork=1:
j-1
sum2=sum2+l(i,k)*l(j,k);
end
l(i,j)=(A(i,j)-sum2)/l(j,j);
end
end
l
y
(1)=b
(1)/l(1,1);
fori=2:
n
sum3=0;
fork=1:
i-1
sum3=sum3+l(i,k)*y(k);
end
y(i)=(b(i)-sum3)/l(i,i);
end
y
x(n)=y(n)/l(n,n);
fori=n-1:
-1:
1
sum4=0;
fork=i+1:
n
sum4=sum4+l(k,i)*x(k);
end
x(i)=(y(i)-sum4)/l(i,i);
end
x
解得的结果如下:
X=[-1.05561.49771.7230-2.71781.9755-0.1730-0.99350.7550-0.57971.1230-0.8631-0.09981.3020-1.29060.1587-0.09921.1349-1.0325-0.11390.8589-0.7735-0.11680.2456-0.05510.9460-2.08892.0415-0.7527-2.40612.3776]’
改进平方根解法程序如下:
n=size(A,1);
d
(1)=A(1,1);
fori=2:
n
forj=1:
i-1
sum1=0;
fork=1:
j-1
sum1=sum1+t(i,k)*l(j,k);
end
t(i,j)=A(i,j)-sum1;
l(i,j)=t(i,j)/d(j);
end
sum2=0;
fork=1:
i-1
sum2=sum2+t(i,k)*l(i,k);
end
d(i)=A(i,i)-sum2;
end
fori=1:
n
l(i,i)=1;
end
y1
(1)=b
(1);
fori=2:
n
sum3=0;
fork=1:
i-1
sum3=sum3+l(i,k)*y1(k);
end
y1(i)=b(i)-sum3;
end
y1
x1(n)=y1(n)/d(n);
fori=n-1:
-1:
1
sum4=0;
fork=i+1:
n
sum4=sum4+l(k,i)*x1(k);
end
x1(i)=y1(i)/d(i)-sum4;
end
x1
解得结果为:
X1=[-1.05561.49771.7230-2.71781.9755-0.1730-0.99350.7550-0.57971.1230-0.8631-0.09981.3020-1.29060.1587-0.09921.1349-1.0325-0.11390.8589-0.7735-0.11680.2456-0.05510.9460-2.08892.0415-0.7527-2.40612.3776]’
实验2.3病态问题的求解
实验目的病态方程的求解更为困难,采用合适的预处理技术可改进问题的病态性
实验题以著名的病态矩阵Hilbert矩阵为例,见本章例6.2,其条件数随阶数
的增加迅速增加,其逆矩阵
,
。
实验步骤和要求
(1)、以下程序作出了
关于阶数
之间的曲线,曲线如下:
clear;
n=30;
%产生hilbert矩阵
fori=1:
1:
n
forj=i:
1:
(n+i-1)
h(i,j-i+1)=1/j;
end
end
forn=1:
1:
100
fori=1:
1:
n
forj=i:
1:
(n+i-1)
h(i,j-i+1)=1/j;
end
end
x(n)=cond(h);
y(n)=log(x(n));
end
%画出曲线
plot(y);
画出的曲线如下如:
(2)画
关于阶数
的曲线的程序如下:
clear;
n=30;
%产生hilbert矩阵
fori=1:
1:
n
forj=i:
1:
(n+i-1)
h(i,j-i+1)=1/j;
end
end
fori=1:
1:
n
D(i,i)=h(i,i);
end
hn=inv(D)*h*inv(D);
x(n)=cond(hn);
y(n)=log(x(n));
plot(y);
图形如下:
(3)、画出
关于阶数
的曲线的程序如下:
clear;
n=30;
%产生hilbert矩阵
fori=1:
1:
n
forj=i:
1:
(n+i-1)
h(i,j-i+1)=1/j;
end
end
fori=1:
1:
n
D(i,i)=h(i,i);
end
hn=inv(D)*h*inv(D);
x(n)=cond(hn)/cond(h);
y(n)=log(x(n));
plot(y);
图形如下:
(4)用高斯消元法解
的程序如下:
clear;
n=6;
%产生hilbert矩阵
fori=1:
1:
n
forj=i:
1:
(n+i-1)
h(i,j-i+1)=1/j;
end
end
b1=h(:
1);
A=h;
b=b1;
na=size(A,2);
nb=size(b,1);
C=[A,b];
ifna~=nb
error('AºÍb±ØÐëÓÐÒ»ÖµÄάÊý');
end
%消元过程
fori=1:
1:
na-1
forj=i+1:
1:
na
C(j,:
)=C(j,:
)-C(j,i)/C(i,i)*C(i,:
);
end
end
fori=1:
1:
na
forj=1:
1:
na
ifC(i,j)<10*eps
C(i,j)=0;
end
end
end
A=C(:
1:
na);
b=C(:
na+1);
%回带求解
bb=b;
fori=na:
-1:
1
x1(i)=bb(i)/A(i,i);
ifi>1
forj=i-1:
-1:
1
bb(j)=bb(j)-A(j,i)*x1(i);
end
end
end
fori=1:
1:
6
D(i,i)=h(i,i);
end
h6=inv(D)*h*inv(D);
b1=h(:
1);
A=h6;
b=b1;
na=size(A,2);
nb=size(b,1);
C=[A,b];
ifna~=nb
error('AºÍb±ØÐëÓÐÒ»ÖµÄάÊý');
end
%消元过程
fori=1:
1:
na-1
forj=i+1:
1:
na
C(j,:
)=C(j,:
)-C(j,i)/C(i,i)*C(i,:
);
end
end
fori=1:
1:
na
forj=1:
1:
na
ifC(i,j)<10*eps
C(i,j)=0;
end
end
end
A=C(:
1:
na);
b=C(:
na+1);
%回带求解
bb=b;
fori=na:
-1:
1
x2(i)=bb(i)/A(i,i);
ifi>1
forj=i-1:
-1:
1
bb(j)=bb(j)-A(j,i)*x2(i);
end
end
end
解得的结果为:
X1=[1,0,0,0,0,0]’
X2=[11,-46.6667,134.4,-205.7143,155.5556,-45.8182]
平方根法解方程
的程序如下:
clear;
n=6;
%产生hilbert矩阵
fori=1:
1:
n
forj=i:
1:
(n+i-1)
h(i,j-i+1)=1/j;
end
end
b1=h(:
1);
A=h;
b=b1;
n=size(A,1);
l(1,1)=sqrt(A(1,1));
fori=2:
n
l(i,1)=A(i,1)/l(1,1);
end
forj=2:
n
sum1=0;
fork=1:
j-1
sum1=sum1+l(j,k)*l(j,k);
end
l(j,j)=sqrt(A(j,j)-sum1);
fori=j+1:
n
sum2=0;
fork=1:
j-1
sum2=sum2+l(i,k)*l(j,k);
end
l(i,j)=(A(i,j)-sum2)/l(j,j);
end
end
l
y
(1)=b
(1)/l(1,1);
fori=2:
n
sum3=0;
fork=1:
i-1
sum3=sum3+l(i,k)*y(k);
end
y(i)=(b(i)-sum3)/l(i,i);
end
x1(n)=y(n)/l(n,n);
fori=n-1:
-1:
1
sum4=0;
fork=i+1:
n
sum4=sum4+l(k,i)*x1(k);
end
x1(i)=(y(i)-sum4)/l(i,i);
end
fori=1:
1:
6
D(i,i)=h(i,i);
end
h6=inv(D)*h*inv(D);
b1=h(:
1);
A=h6;
b=b1;
n=size(A,1);
l(1,1)=sqrt(A(1,1));
fori=2:
n
l(i,1)=A(i,1)/l(1,1);
end
forj=2:
n
sum1=0;
fork=1:
j-1
sum1=sum1+l(j,k)*l(j,k);
end
l(j,j)=sqrt(A(j,j)-sum1);
fori=j+1:
n
sum2=0;
fork=1:
j-1
sum2=sum2+l(i,k)*l(j,k);
end
l(i,j)=(A(i,j)-sum2)/l(j,j);
end
end
y
(1)=b
(1)/l(1,1);
fori=2:
n
sum3=0;
fork=1:
i-1
sum3=sum3+l(i,k)*y(k);
end
y(i)=(b(i)-sum3)/l(i,i);
end
x2(n)=y(n)/l(n,n);
fori=n-1:
-1:
1
sum4=0;
fork=i+1:
n
sum4=sum4+l(k,i)*x2(k);
end
x2(i)=(y(i)-sum4)/l(i,i);
end
解得的结果为:
X1=[1,0,0,0,0,0];
X2=[11,-46.6667,134.4,-205.7143,155.5556,-45.8182]
实验2.4主元的选取与算法稳定性
实验目的主元的选择从理论上看很平凡,却是数值计算中十分典型的问题
实验题系数矩阵
,方程组的精确解为
实验步骤与要求
(1)自行编一个能自动选主元,又能手动选取主元的求解线性方程组的Gauss消元过程的程序;
(2)取阶数
,用Matlab的“condest”估计矩阵条件数;(3)用程序自动选主元,计算结果与精确解作比较;
(1)自主选主元的高斯消元法程序如下:
C=[Ab];
n=size(A,2);
fori=2:
n
[r,c]=max(abs(C(i-1:
n,i-1)));
temp=C(i-1,:
);
C(i-1,:
)=C(c+i-2,:
);
C(c+i-2,:
)=temp;
ifC(i-1,i-1)==0
return
end
fork=(i-1):
(n-1)
C(k+1,:
)=C(i-1,:
)*-C(k+1,i-1)/C(i-1,i-1)+C(k+1,:
);
end
end
D=C;
B=C(:
n+1);
D(:
n+1)=[];
x=zeros(1,n);
fori=1:
n
x(n+1-i)=(B(n+1-i)-D(n+1-i,:
)*x')/C(n+1-i,n+1-i);
end
(2)、产生阶数分别为n=10,30,100的矩阵A以及估计各个矩阵条件数的程序如下:
clear;
%n=10的系数矩阵A
n=10;
fori=1:
1:
n
if(i==1)
A10(i,i)=6;
A10(i,i+1)=1;
elseif(i==n)
A10(i,i-1)=8;
A10(i,i)=6;
else
A10(i,i-1)=8;
A10(i,i)=6;
A10(i,i+1)=1;
end
end
end
%n=30的系数矩阵A
n=30;
fori=1:
1:
n
if(i==1)
A30(i,i)=6;
A30(i,i+1)=1;
elseif(i==n)
A30(i,i-1)=8;
A30(i,i)=6;
else
A30(i,i-1)=8;
A30(i,i)=6;
A30(i,i+1)=1;
end
end
end
%n=100的系数矩阵A
n=100;
fori=1:
1:
n
if(i==1)
A100(i,i)=6;
A100(i,i+1)=1;
elseif(i==n)
A100(i,i-1)=8;
A100(i,i)=6;
else
A100(i,i-1)=8;
A100(i,i)=6;
A100(i,i+1)=1;
end
end
end
m1=condest(A10);
m2=condest(A30);
m3=condest(A100);
得结果如下:
(3)、对于n=10的矩阵,由
(1)中的自主选主法得其结果为
X10=[1,1,1,1,1,1,1,1,1],与精确解一致
对于n=30的矩阵,由
(1)中的自主选主法所得的结果任然和精确解一致
对于n=100的矩阵,由
(1)中的自主选主法所得的结果如下:
x