数值分析实验报告二
求解线性方程组的直接方法(2学时)
班级专业 统计一班 姓名 何斌 学号 06 日期 2011-3-30
一 实验目的
1.掌握求解线性方程组的高斯消元法及列主元素法; 2. 掌握求解线性方程组的克劳特法; 3. 掌握求解线性方程组的平方根法。
二 实验内容
1.用高斯消元法求解方程组(精度要求为??10?6):
?3x1?x2?2x3?7???x1?2x2?2x3??1 ?2x?2x?4x?023?12.用克劳特法求解上述方程组(精度要求为??10?6)。 3. 用平方根法求解上述方程组(精度要求为??10?6)。 4. 用列主元素法求解方程组(精度要求为??10?6):
?3x1?x2?4x3?3???x1?2x2?2x3?2 ?2x?3x?2x??523?1三 实验步骤(算法)与结果
1. 高斯消元法求解:
设所求方程组的增广矩阵为A0,系数矩阵为A01,(择定值)L=(l11,l22,l33), 方程组右边的值用有序数组B表示. clear
A0=[3 -1 2 7;-1 2 -2 -1;2 -2 4 0];L=[1,1,1]; A01=[3 -1 2 ;-1 2 -2 ;2 -2 4 ];B=[7 -1 0]; %a=zeros(size(A0)); a=A0;
%l=zeros(1,3); l=L;
1
%消元过程的计算公式
u11=a(1,1)/l(1,1);u12=a(1,2)/l(1,1);u13=a(1,3)/l(1,1);z1=a(1,4)/l(1,1); a(2,1)=a(2,1)/u11;a(2,2)=(a(2,2)-a(2,1)*u12)/l(1,2);
a(2,3)=(a(2,3)-a(2,1)*u13)/l(1,2);z2=(a(2,4)-a(2,1)*z1)/l(1,2); a(3,1)=a(3,1)/u11;a(3,2)=(a(3,2)-a(3,1)*u12)/a(2,2); a(3,3)=(a(3,3)-a(3,1)*u13-a(3,2)*a(2,3))/l(1,3); z3=(a(3,4)-a(3,1)*z1-a(3,2)*z2)/l(1,3);
A1=[u11 u12 u13 ;0 a(2,2) a(2,3) ; 0 0 a(3,3)]; %得到上三角矩阵 Z=[z1;z2;z3];
%迭代求解
x3=z3/a(3,3); x2=(z2-a(2,3)*x3)/a(2,2); x1=(z1-u13*x3-u12*x2)/u11; x=[x1 x2 x3]; for k=1:3
fprintf('x(%d)的结果为%.6f\\n\\n',k,x(k)),end %控制精度 %残差法检验答案 C=A01*X'-B'
程序运行结果:
x(1)的结果为3.500000 x(2)的结果为-1.000000 x(3)的结果为-2.250000 残差法检验结果:
C = 1.0e-014 *[ 0.0888, -0.0888, 0.1776]’ 显然符合要求
2.克劳特法求解:
易知克劳特法与高斯消元法求解中仅择取的L初值不同,故改取L=[3 2 4]即可。运行上述程序得: X1=3.500000 x2=-1. x3=-2.250000 C =[0 0 0] 显然符合要求
3.平方根法求解:
易知改取L=[3^(1/2) 2^(1/2) 4^(1/2)]即可。运行上述程序得: X1=3.500000 x2=-1. x3=-2.250000
C = 1.0e-014 *[ 0.0888, 0, 0.1776]’ 显然符合要求
4. 列主元素法: clear
%列主元素法
A0=[3 -1 4 3;-1 2 -2 2;2 -3 -2 -5];%原方程组的增广矩阵
[m,n]=max(A0(:,1));%查找增广矩阵A0中第一列中最大的数m,得其行号为n if n~=1 %最大数不在第一行,则进行行的转换 b1=A0(n,:);%提取原方程第n行x的系数 b11=zeros(size(A0)); b11(1,:)=b1;
a1=A0(1,:);% 提取原方程第一行x的系数
2
a11=zeros(size(A0)); a11(n,:)=a1;
A0(n,:)=[0 0 0 0]; A0(1,:)=[0 0 0 0];
A0=A0+b11+a11;%实现了行的转换,即第一行与第n行 end a=A0;
%计算乘数
l21=a(2,1)/a(1,1);l31=a(3,1)/a(1,1);
%消元
a(2,:)=a(2,:)-l21*a(1,:);a(3,:)=a(3,:)-l31*a(1,:); aa=a([2 3],[2 3 4]); [m,n]=max(aa(:,1));
if n~=1 circshift(aa,1);end %最大数不在第一行,则进行行的转换 %计算乘数
l32=aa(2,1)/aa(1,1);
%消元
aa(2,:)=aa(2,:)-l32*aa(1,:);
%回代求解
x3=aa(2,3)/aa(2,2);x2=(aa(1,3)-aa(1,2)*x3)/aa(1,1); x1=(a(1,4)-a(1,3)*x3-a(1,2)*x2)/a(1,1);x=[x1 x2 x3]; for k=1:3
fprintf('x(%d)的结果为%.6f\\n\\n',k,x(k)),end %残差法检验
A01=A0([1 2 3],[1 2 3]); B=A0(:,4);
C=A01*[x1 x2 x3]'-B
程序运行结果:
x(1)的结果为1.000000 x(2)的结果为2.000000 x(3)的结果为0.500000 C =[0 0 0]’ 显然符合要求
四 实验收获与教师评语
这次实验使我对MATLAB的数组运算有了更深的了解,程序得到了很好的实现,能方便地解三元一次线性方程组,但是本程序仍然有很大的局限性。
3

