数值计算方法 课后实验 李维国 中国石油大学
合集下载
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
m=A(i,k)/A(k,k);
for j=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
if abs(A(n,n))<1e-10
index=0;return;
end
for k=n:-1:1
（4）观察小主元并分析对计算结果的影响．
解：（1）
程序如下：
clear
>> A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
>> cond(A)
结果如下：
ans =
68.4296
分析：得到的结果cond（A）=68.4296>10所以该矩阵属于病态。
end
if k==n
break;
end
lu(k+1:n,k)=lu(k+1:n,k)/lu(k,k);
lu(k+1:n,k+1:n)=lu(k+1:n,k+1:n)-lu(k+1:n,k)*lu(k,k+1:n);
end
l=diag(ones(n,1))+tril(lu,-1);
u=triu(lu);
y(1)=b(p(1));
for i=2:n
y(i)=b(p(i))-l(i,1:i-1)*y(1:i-1)';
end
x(n)=y(n)/u(n,n);
for i=n-1:-1:1
x(i)=(y(i)-u(i,i+1:n)*x(i+1:n)')/u(i,i);
end
x=x';
执行程序：clear
A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
function x=gauss_lie(A,b)
%采用高斯列主元法求解方程组Ax=b
n=length(b);
p=1:n;lu=A;
y=[];
for k=1:n
[c,i]=max(abs(lu(k:n,k)));
ik=i+k-1;
if ik~=k
m=p(k);p(k)=p(ik);p(ik)=m;
ck=lu(k,:);lu(k,:)=lu(ik,:);lu(ik,:)=ck;
>> [L,U,P]=lu(A)
L =
1.0000 0 0 0
0.0000 1.0000 0 0
0.4724 -0.1755 1.0000 0
0.0893 0.0202 -0.1738 1.0000
U =
11.2000 9.0000 5.0000 2.0000
0 59.1400 3.0000 1.0000
for i=k:n
if abs(A(i,k))a_max
a_max=abs(A(i,k));r=i;
end
end
if r>k
for j=k:n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=bຫໍສະໝຸດ Baiduk);b(k)=b(r);b(r)=z;det=-det;
end
for i=k+1:n
det =
-1.9065e+003
index =
1
（3）不选主元的分解程序如下：
clear
A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
>> [l,u]=lu(A)
得到LU分解结果如下：
l =
0.0000 1.0000 0 0
b=[59.17,46.78,1,2];
>> x=gauss_lie(A,b)
结果如下：得到的x向量
x =
3.8457
1.6095
-15.4761
10.4113
（3）选列主元的分解如下：
clear
>> A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
0.4724 -0.1755 1.0000 0
1.0000 0 0 0
0.0893 0.0202 -0.1738 1.0000
u =
11.2000 9.0000 5.0000 2.0000
0 59.1400 3.0000 1.0000
0 0 -2.8354 1.2307
0 0 0 1.0151
（2）高斯列主元消去法二：
实验3.1 Gauss消去法的数值稳定性实验
实验目的：观察的理解高斯消元过程中出现小主元即很小时，引起方程组解数值不稳定性．
实验题：求解线性方程组．
，
实验要求：
（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的．
（2）用高斯列主元消去法求得L和U及解向量．
（3）用不选主元的高斯消去法求得L和U及解向量．
（2）高斯列主元消去法程序如下：
function [x,det,index]=gauss3(A,b)
%x为方程组的解，det为系数行列式的值，index标志位，值为1表示执行成功，否则执行失败。
n=length(b);
index=1;det=1;
x=zeros(1,n);
for k=1:n-1
a_max=0;
for j=k+1:n
b(k)=b(k)-A(k,j)*x(j);
end
x(k)=b(k)/A(k,k);
end
执行下列程序：
clear
A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
b=[59.17,46.78,1,2];
[x,det,index]=gauss3(A,b)
得到结果如下：
a_max =
0
a_max =
3.0000e-016
a_max =
5.2910
a_max =
11.2000
a_max =
0
a_max =
16.7120
a_max =
13.4000
a_max =
0
a_max =
5.5203
x =
3.8457 1.6095 -15.4761 10.4113
0 0 -2.8354 1.2307
0 0 0 1.0151
P =
0 0 1 0
1 0 0 0
0 1 0 0
0 0 0 1
for j=k+1:n
A(i,j)=A(i,j)-m*A(k,j);
end
b(i)=b(i)-m*b(k);
end
det=det*A(k,k);
end
det=det*A(n,n);
if abs(A(n,n))<1e-10
index=0;return;
end
for k=n:-1:1
（4）观察小主元并分析对计算结果的影响．
解：（1）
程序如下：
clear
>> A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
>> cond(A)
结果如下：
ans =
68.4296
分析：得到的结果cond（A）=68.4296>10所以该矩阵属于病态。
end
if k==n
break;
end
lu(k+1:n,k)=lu(k+1:n,k)/lu(k,k);
lu(k+1:n,k+1:n)=lu(k+1:n,k+1:n)-lu(k+1:n,k)*lu(k,k+1:n);
end
l=diag(ones(n,1))+tril(lu,-1);
u=triu(lu);
y(1)=b(p(1));
for i=2:n
y(i)=b(p(i))-l(i,1:i-1)*y(1:i-1)';
end
x(n)=y(n)/u(n,n);
for i=n-1:-1:1
x(i)=(y(i)-u(i,i+1:n)*x(i+1:n)')/u(i,i);
end
x=x';
执行程序：clear
A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
function x=gauss_lie(A,b)
%采用高斯列主元法求解方程组Ax=b
n=length(b);
p=1:n;lu=A;
y=[];
for k=1:n
[c,i]=max(abs(lu(k:n,k)));
ik=i+k-1;
if ik~=k
m=p(k);p(k)=p(ik);p(ik)=m;
ck=lu(k,:);lu(k,:)=lu(ik,:);lu(ik,:)=ck;
>> [L,U,P]=lu(A)
L =
1.0000 0 0 0
0.0000 1.0000 0 0
0.4724 -0.1755 1.0000 0
0.0893 0.0202 -0.1738 1.0000
U =
11.2000 9.0000 5.0000 2.0000
0 59.1400 3.0000 1.0000
for i=k:n
if abs(A(i,k))a_max
a_max=abs(A(i,k));r=i;
end
end
if r>k
for j=k:n
z=A(k,j);A(k,j)=A(r,j);A(r,j)=z;
end
z=bຫໍສະໝຸດ Baiduk);b(k)=b(r);b(r)=z;det=-det;
end
for i=k+1:n
det =
-1.9065e+003
index =
1
（3）不选主元的分解程序如下：
clear
A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
>> [l,u]=lu(A)
得到LU分解结果如下：
l =
0.0000 1.0000 0 0
b=[59.17,46.78,1,2];
>> x=gauss_lie(A,b)
结果如下：得到的x向量
x =
3.8457
1.6095
-15.4761
10.4113
（3）选列主元的分解如下：
clear
>> A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
0.4724 -0.1755 1.0000 0
1.0000 0 0 0
0.0893 0.0202 -0.1738 1.0000
u =
11.2000 9.0000 5.0000 2.0000
0 59.1400 3.0000 1.0000
0 0 -2.8354 1.2307
0 0 0 1.0151
（2）高斯列主元消去法二：
实验3.1 Gauss消去法的数值稳定性实验
实验目的：观察的理解高斯消元过程中出现小主元即很小时，引起方程组解数值不稳定性．
实验题：求解线性方程组．
，
实验要求：
（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的．
（2）用高斯列主元消去法求得L和U及解向量．
（3）用不选主元的高斯消去法求得L和U及解向量．
（2）高斯列主元消去法程序如下：
function [x,det,index]=gauss3(A,b)
%x为方程组的解，det为系数行列式的值，index标志位，值为1表示执行成功，否则执行失败。
n=length(b);
index=1;det=1;
x=zeros(1,n);
for k=1:n-1
a_max=0;
for j=k+1:n
b(k)=b(k)-A(k,j)*x(j);
end
x(k)=b(k)/A(k,k);
end
执行下列程序：
clear
A=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1];
b=[59.17,46.78,1,2];
[x,det,index]=gauss3(A,b)
得到结果如下：
a_max =
0
a_max =
3.0000e-016
a_max =
5.2910
a_max =
11.2000
a_max =
0
a_max =
16.7120
a_max =
13.4000
a_max =
0
a_max =
5.5203
x =
3.8457 1.6095 -15.4761 10.4113
0 0 -2.8354 1.2307
0 0 0 1.0151
P =
0 0 1 0
1 0 0 0
0 1 0 0
0 0 0 1