数值分析编程及运行结果(高斯顺序消元法)

数值分析实验,用程序实现Gauss消元法

《数值分析》实验报告实验序号：实验一实验名称:一实验目的：用Microsoft visual c++编写一个功能与Windows 计算器一样的计算器。

二实验内容：(1)算法介绍：计算方法分析：Gauss消元法的基本做法就是把方程组转化成为一个如下图的等价的三角方程组，这个过程叫做消元。

得到三角方程组后，就可以逐个求出Xn，Xn-1，…，X1，这个过程叫回代。

程序代码分析：建立两个数组a和b，通过循环语句将ｎ阶增广矩阵输入进去，通过对列的循环对每一列进行消去未知数，通过ｎ小步ｎ大步把矩阵化简成上三角形矩阵，最后通过迭代法解得方程组得解。

(2)算法分析：程序中只用到一个主函数，求解线形方程组得算法都放在主函数中，利用以下函数进行求解：a[i][j] = a[i][j] - (a[i][k] / a[k][k]) * a[k][j];迭代：x[n] = a[n][n + 1] / a[n][n];for(i = n - 1;i >= 1;i --){sum = 0.0;for(j = n;j >= i + 1;j --){sum = sum + a[i][j] * x[j];}//cout << sum << endl;x[i] = (a[i][n + 1] - sum) / a[i][i];}(3)原代码：#include<iostream>#include <stdio.h>#include <cmath>#define N 1000using namespace std;double a[N][N];double b[N];int main(){int n,i,j,k;double m;double x[N];printf("***Gauss消元法***:\n\n");while(printf("输入矩阵的阶数：")){scanf("%d",&n);printf("输入增广矩阵：\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++){scanf("%lf",&a[i][j]); ;}scanf("%lf",&b[i]);}for(i=1;i<n;i++){for(k=i+1;k<=n;k++){m=a[k][i]*1/a[i][i]*1;for(j=i;j<=n;j++){a[k][j]= a[k][j]-a[i][j]*m;}b[k]= b[k]-b[i]*m*1;}printf("经过第%d大步后增广矩阵为：\n",i);for(int l=1;l<=n;l++){for(j=1;j<=n;j++)printf("%f ",a[l][j]);printf("%f ",b[l]);printf("\n");}}for(i=n;i>0;i--){double sum=0;for(j=i+1;j<=n;j++){sum =sum+a[i][j]*x[j]*1;}x[i]=(b[i]-sum)*1/a[i][i]*1;}printf("线性方程组的解为: \n");for(i=1;i<=n;i++){if(x[i]==0)x[i]=abs(x[i]);printf("x[%d]=%f\n",i,x[i]);}}return 0;}三实验截图：四实验结果分析：消去第k个元素时，对矩阵作加法和乘法各（n-k）*（n-k）次，除法运算（n-k）次，对右端作加法和乘法各（n-k）次，加法和乘法运算分别各n*(n-1)*(n-1/2)/3和n*(n-1)/2次，消元时还有n*(n-1)/2次除法运算，另外回代过程加法和乘法运算各n*(n-1)/2，除法运算n次。

数值分析编程及运行结果(高斯顺序消元法)课案

高斯消元法1.程序:clearformat ratA=input('输入增广矩阵A=')[m,n]=size(A);for i=1:(m-1)numb=int2str(i);disp(['第',numb,'次消元后的增广矩阵']) for j=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endAend%回代过程disp('回代求解')x(m)=A(m,n)/A(m,m);for i=(m-1):-1:1x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); endx2.运行结果：高斯选列主元消元法1.程序:clearformat ratA=input('输入增广矩阵A=')[m,n]=size(A);for i=1:(m-1)numb=int2str(i);disp(['第',numb,'次选列主元后的增广矩阵']) temp=max(abs(A(i:m,i)));[a,b]=find(abs(A(i:m,i))==temp);tempo=A(a(1)+i-1,:);A(a(1)+i-1,:)=A(i,:);A(i,:)=tempodisp(['第',numb,'次消元后的增广矩阵'])for j=(i+1):mA(j,:)=A(j,:)-A(i,:)*A(j,i)/A(i,i);endAend%回代过程disp('回代求解')x(m)=A(m,n)/A(m,m);for i=(m-1):-1:1x(i)=(A(i,n)-A(i,i+1:m)*x(i+1:m)')/A(i,i); endx2.运行结果：追赶法1.程序：function [x,L,U]=zhuiganfa(a,b,c,f)a=input('输入矩阵-1对角元素a=');b=input('输入矩阵对角元素b=');c=input('输入矩阵+1对角元素c=');f=input('输入增广矩阵最后一列元素f='); n=length(b);% 对A进行分解u(1)=b(1);for i=2:nif(u(i-1)~=0)l(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);elsebreak;endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;% 求解Ly=by(1)=f(1);for i=2:ny(i)=f(i)-l(i-1)*y(i-1); end% 求解Ux=yif(u(n)~=0)x(n)=y(n)/u(n);endfor i=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1))/u(i); end2.运行结果：高斯-塞德尔迭代格式1.程序：function x=Gauss_Seidel(a,b)a=input('输入系数矩阵a=')b=input('输入增广矩阵最后一列b=');e=0.5e-7;n=length(b);N=50;x=zeros(n,1);t=zeros(n,1);for k=1:Nsum=0;E=0;t(1:n)=x(1:n);for i=1:nx(i)=(b(i)-a(i,1:(i-1))*x(1:(i-1))-a(i,(i+1):n)*t((i+1):n))/a(i,i);endif norm(x-t)<ekbreak;endend2.运行结果：雅戈比迭代格式1.程序：function x=Jocabi(a,b)a=input('输入系数矩阵a=');b=input('输入增广矩阵最后一列b=');e=0.5e-7;n=length(b);N=100;x=zeros(n,1);y=zeros(n,1);for k=1:Nsum=0;for i=1:ny(i)=(b(i)-a(i,1:n)*x(1:n)+a(i,i)*x(i))/a(i,i);endfor i=1:nsum=sum+(y(i)-x(i))^2;endif sqrt(sum)<ekbreak;elsefor i=1:nx(i)=y(i);endendendif k==N warning('未能找到近似解'); end2.运行结果：逐次超松弛法（SOR）1.程序：function [n,x]=sor22(A,b,X,nm,w,ww)%用超松弛迭代法求解方程组Ax=b%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子%输出：x为求得的方程组的解构成的列向量，n为迭代次数A=input('输入系数矩阵A=');b=input('输入方程组右端的列向量b=');X=input('输入迭代初值构成的列向量X=');nm=input('输入最大迭代次数nm=');w=input('输入误差精度w=');ww=input('输入松弛因子ww=');n=1;m=length(A);D=diag(diag(A)); %令A=D-L-U,计算矩阵DL=tril(-A)+D; %令A=D-L-U,计算矩阵LU=triu(-A)+D; %令A=D-L-U,计算矩阵UM=inv(D-ww*L)*((1-ww)*D+ww*U); %计算迭代矩阵g=ww*inv(D-ww*L)*b; %计算迭代格式中的常数项%下面是迭代过程while n<=nmx=M*X+g; %用迭代格式进行迭代if norm(x-X,'inf')<wdisp('迭代次数为');ndisp('方程组的解为');xreturn;%上面：达到精度要求就结束程序，输出迭代次数和方程组的解endX=x;n=n+1;end%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）disp('在最大迭代次数内不收敛!');disp('最大迭代次数后的结果为');2.运行结果：二分法求解方程的根1.程序：%其中a，b表示查找根存在的范围，M表示要求解函数的值%f(x)表示要求解根的方程%eps表示所允许的误差大小function y=er_fen_fa(a,b,M)k=0;eps=0.05while b-a>epsx=(a+b)/2;%检查是否大于值if ((x^3)-3*x-1)>Mb=xelsea=xendk=k+1end2.运行结果：Newton 迭代法（切线法）1.程序：function x=nanewton(fname,dfname,x0,e,N)%newton迭代法解方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内嵌函数，x0为迭代初值，e为精度要求x=x0;x0=x+2*e;k=0;if nargin<5,N=500;endif nargin<4 e=1e-4;endwhile abs(x0-x)>e&k<N,k=k+1;x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);disp(x)endif k==N,warning('已达迭代次数上限');end2.运行结果：割线方式迭代法1.程序：function x=ge_xian_fa(fname,dfname,x0,x1,e,N)%割线方式迭代法解方程组%fname和dfname分别表示F(x)及其导函数的M函数句柄或内嵌函数，x0,x1分别为迭代初值，e为精度要求k=0;a=x1;b=x0;if nargin<5,N=500;endif nargin<4 e=1e-4;endwhile abs(a-b)>e&k<N,k=k+1;x=x1-((x1-x0)/(feval(fname,x1)-feval(fname,x0)))*feval(fname,x1);if feval(fname,x)*feval(fname,x0)>0,x0=x;b=x0;elsex1=x;a=x1;endx=x1-((x1-x0)/(feval(fname,x1)-feval(fname,x0)))*feval(fname,x1);numb=int2str(k);disp(['第',numb,'次计算后x='])fprintf('%f\n\n',x);endif k==N,warning('已达迭代次数上限'); end2.运行结果：Newton插值1.程序：%保存文件名为New_Int.m%Newton基本插值公式%x为向量，全部的插值节点%y为向量，差值节点处的函数值%xi为标量，是自变量%yi为xi出的函数估计值function yi=newton_chazhi(x,y,xi)n=length(x);m=length(y);if n~=merror('The lengths of X ang Y must be equal!'); return;end%计算均差表YY=zeros(n);Y(:,1)=y';for k=1:n-1for i=1:n-kif abs(x(i+k)-x(i))<epserror('the DATA is error!');return;endY(i,k+1)=(Y(i+1,k)-Y(i,k))/(x(i+k)-x(i)); endend%计算牛顿插值公式yi=0;for i=1:nz=1;for k=1:i-1z=z*(xi-x(k));endyi=yi+Y(1,i)*z;end2.运行结果：Lagrange插值1.程序：function y0 = Language(x,y,x0)syms t l;if length(x)==length(y)n = length(x);elsedisp('x和y的维数不相等！');return; %检错endh=sym(0);for i=1:nl=sym(y(i));for j=1:i-1l=l*(t-x(j))/(x(i)-x(j));end;for j=i+1:nl=l*(t-x(j))/(x(i)-x(j));end;h=h+l;endsimplify(h);if nargin == 3y0 = subs (h,'t',x0); %计算插值点的函数值elsey0=collect(h);y0 = vpa(y0,6); %将插值多项式的系数化成6位精度的小数end2.运行结果：最小二乘法1.程序：function p=nafit(x,y,m)%多项式拟合% x, y为已知数据点向量, 分别表示横,纵坐标, m为拟合多项式的次数, 结果返回m次拟合多项式系数, 从高次到低次存放在向量p 中.A=zeros(m+1,m+1);for i=0:mfor j=0:mA(i+1,j+1)=sum(x.^(i+j));endb(i+1)=sum(x.^i.*y);enda=A\b';p=fliplr(a');t=0:0.1:1.6;S3=polyval(p,t);plot(x,y,'p k');hold onplot(t,S3,'r');2.运行结果：。

数值分析计算方法实验报告

break;
end;
end;
X=x;
disp('迭代结果:');
X
format short;
输出结果：
因为不收敛，故出现上述情况。
4.超松弛迭代法：
%SOR法求解实验1
%w=1.45
%方程组系数矩阵
clc;
A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]
b=[10,5,-2,7]'
b=[10,5,-2,7]'
[m,n]=size(A);
if m~=n
error('矩阵A的行数和列数必须相同');
return;
end
if m~=size(b)
error('b的大小必须和A的行数或A的列数相同');
return;
end
if rank(A)~=rank([A,b])
error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');
3．实验环境及实验文件存档名
写出实验环境及实验文件存档名
4.实验结果及分析
输出计算结果，结果分析和小结等。
解：1.高斯列主元消去法：
%用高斯列主元消去法解实验1
%高斯列主元消元法求解线性方程组Ax=b
%A为输入矩阵系数,b为方程组右端系数
%方程组的解保存在x变量中
format long;
A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]
return;
end
c=n+1;
A(:,c)=b;
for k=1:n-1

数值分析列主元高斯消去法

实验四：列组元消去法一、目的1）熟悉列主元高斯消元法解线性方程组的算法2）掌握列主元高斯消去法的编程二、实验原理列主元素消去法是为控制舍入误差而提出来的一种算法,在Gauss消去法的消元过程中,若出现a=0,则消元无法进行,即使其不为0,但很小,把它作为除数,就会导致其他元素量级的巨大增长和舍入误差的扩散,最后使计算结果不可靠.使用列主元素消去法计算,基本上能控制舍入误差的影响,并且选主元素比较方便.三、运行结果四、代码using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace高斯{class Program{static double[] Gause(double[,] a, int n){int i, j, k;int rank, columm;double temp, l, s, mx;double[] x = new double[n];for (i = 0; i <= n - 2; i++){mx = Math.Abs(a[i, i]);rank = i;columm = i;for (j = i + 1; j <= n - 1; j++) //选主元if (Math.Abs(a[j, i]) > mx){mx = Math.Abs(a[j, i]);rank = j;columm = i;}for (k = 0; k <= n; k++) //主元行变换{temp = a[i, k];a[i, k] = a[rank, k];a[rank, k] = temp;} //消元for (j = i + 1; j <= n - 1; j++){l = a[j, i] / a[i, i];for (k = i; k <= n; k++)a[j, k] = a[j, k] - l * a[i, k];}}x[n - 1] = a[n - 1, n] / a[n - 1, n - 1]; //回代方程求解x for (i = n - 2; i >= 0; i--){s = 0;for (j = i + 1; j <= n - 1; j++)s = s + a[i, j] * x[j];x[i] = (a[i, n] - s) / a[i, i];}return x;}static void Main(string[] args){double[,] a = new double[4, 5] { { 10, -7, 0, 1, 8 }, { -3, 2.099999, 6, 2, 5.900001 }, { 5, -1, 5, -1, 5 }, { 2, 1, 0, 2, 1 } };int n = 4;double[] x = new double[n];x = Gause(a, n);Console.WriteLine("高斯消去法方程：");for (int i = 0; i < n; i++){for (int j = 0; j < n; j++)Console.Write(a[i, j].ToString() + " ");Console.WriteLine();}Console.WriteLine("线性方程组的解：");for (int i = 0; i <= n - 1; i++)Console.Write("x" + (i + 1).ToString() + "=" +x[i].ToString() + " ");Console.WriteLine();Console.ReadLine();}}}四、分析通过本次实验的学习，学会根据算法编写基本的相关程序，虽然此次程序模板由老师给予，但认真阅读理解程序有助于今后的学习，再利用计算机中的C语言对高斯列主元消去法可以快速得到线性方程组的解，由简单的线性方程组可以推广到一般n阶线性方程组，这对如何利用高斯列主元消去法解决实际问题有了一定的经验。

数值分析4 高斯主元素消去法

§2高斯主元素消去法⎪⎩⎪⎨⎧=++-=++=++00.357.404.100.200.224.563.200.100.100.200.10120.0321321321x x x x x x x x x 解：clear alla=[0.0120 1.00 2.00;1.00 2.63 5.24;-2.00 1.04 4.57]; b=[1.00;2.00;3.00];x=a\b方程组的三位有效数字的解：Tx )266.0,476.0,645.0(*-=Gauss 消去法求解(取三位有效数字)：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---−−→−⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡------−−−→−⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-==-==00.300.5003.811627.80000.100.200.10120.016432916603.811627.80000.100.200.10120.000.357.404.100.200.224.563.200.100.100.200.10120.006.21673.83323121l l l b A 解出Tx )60.0,197.0,25.0(--≈。

【注】1）设Ax=b,其中A 为n 阶非奇异矩阵，可以应用高斯消元法。

2）消元过程中，即使0)(≠k kk a ，用其作除数)/()()(k kk k ik ik a a l =会导致计算中间结果数量级严重增长和舍入误差的累积、扩大，最后使得计算结果不可靠。

3)应避免采用绝对值很小的主元素)(k kk a ；对一般的系数矩阵，最好保持乘数1≤ik l ，因此，在高斯消去法中应引进选主元技巧，以便减少计算过程中舍入误差对求解的影响。

clear alla=[0.0120 1.00 2.00;1.00 2.63 5.24;-2.00 1.04 4.57]; b=[1.00;2.00;3.00];x_value=vpa(a\b,15)%10位有效数字的近似解a=[a,b];eps=1e-6;[n,m]=size(a);Gauss,x=vpa(x,15) %对比高斯消去法的结果一、列主元素消去法基本思想：在每轮消元之前，选列主元素（绝对值最大的元素）,使乘数（即消元因子）1≤ik l步骤：设已进行k-1轮消元，得矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=)()()()()2(2)2(2)2(22)1(1)1(1)1(12)1(11)(k nn k nkk kn k kkn kn kk a a a a a a a a a a a A一落千丈 1 23Ｓ１：选列主元素： )()(0max k ik ni k k k i a a ≤≤= (1)Ｓ２：换行：如果)(0k k i a →０，则方程组解不唯一，停止运算；否则，如果i0=k ，则可进行下一轮消元；如果k i ≠0，则r i0 r k ，然后进行下一轮消元。

数值分析高斯顺序消去法、列主元消去法LU分解法

数值分析实验报告（1）学院：信息学院班级：计算机0903班姓名：***学号：********课题一A.问题提出给定下列几个不同类型的线性方程组，请用适当的方法求解线性方程组1、设线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------------------1368243810041202913726422123417911101610352431205362177586832337616244911315120130123122400105635680000121324⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125 x *= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 )T2、设对称正定阵系数阵线方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------------19243360021411035204111443343104221812334161206538114140231212200420424⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡87654321x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---4515229232060 x * = ( 1, -1, 0, 2, 1, -1, 0, 2 )T3、三对角形线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------4100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----5541412621357 x *= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 )TB.（1）对上述三个方程组分别用Gauss 顺序消去法与Gauss 列主元消去法；平方根与改进平方根法；追赶法求解（选择其一）（2）编写算法通用程序（3）在应用Gauss 消去时，尽可能利用相应程序输出系数矩阵的三角分解式C.(1)通过该课题的程序编制，掌握模块化结构程序设计方法 (2)掌握求解各类线性方程组的直接方法，了解各种方法的特点 (3)体会高斯消去法选主元的必要性实验步骤：（高斯消去法，列主元，LU ）1顺序高斯消去法2.LU 分解法3.列主元高斯消去法（如下图）（1）高斯消去法运行结果如下（2）对方程的系数矩阵进行LU分解并求出方程组的解（3）列主元高斯消去法实验体会总结：利用gauss消去法解线性方程组的时候，如果没有经过选主元，可能会出现数值不稳定的现象，使得方程组的解偏离精确解。

数值分析列主元高斯消去顺序高斯平方根法追赶法

课题名称：课题一解线性方程组的直接方法解决的问题：给定三个不同类型的线性方程组，用适当的直接法求解。

采用的数值方法：对第一个普通的线性方程组，采用了高斯顺序消去法和高斯列主元消去法。

对第二个正定线性方程组，采用了平方根法。

对第三个三对角线性方程组，采用了追赶法。

算法程序：(1)普通的线性方程组①顺序消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0},{8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;int size=10;for(k=0; k<size-1; k++){if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k];for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j];}x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}②列主元消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0}, {8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;float max;int col;int size=10;for(k=0; k<size-1; k++){max=fabs(A[k][k]);col=k;for(i=k; i<size; i++){if(max<fabs(A[i][k])) {max=fabs(A[i][k]); col=i;}}for(j=k; j<size; j++){temp=A[col][j];A[col][j]=A[k][j];A[k][j]=temp;}temp=b[col];b[col]=b[k];b[k]=temp;if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k]; for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j]; }x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}(2)对称正定线性方程组平方根法：#include <stdio.h>#include <math.h>#define n 8int main(void){float A[8][8]={{4,2,-4,0,2,4,0,0},{2,2,-1,-2,1,3,2,0},{-4,-1,14,1,-8,-3,5,6},{0,-2,1,6,-1,-4,-3,3},{2,1,-8,-1,22,4,-10,-3},{4,3,-3,-4,4,11,1,-4},{0,2,5,-3,-10,1,14,2},{0,0,6,3,-3,-4,2,19}};float g[8][8]= {0};float b[8]= {0,-6,6,23,11,-22,-15,45}; float x[8]= {0};float y[8]= {0};int k,m,i,sq;for(k=0; k<n; k++){float p=0,q=0,s=0;for(m=0; m<=k-1; m++){p=p+A[k][m]*A[k][m];}g[k][k]=sqrt(A[k][k]-p);A[k][k]=g[k][k];for(i=k+1; i<n; i++){q=0;for(m=0; m<=k-1; m++){q=q+A[i][m]*A[k][m];}g[i][k]=(A[i][k]-q)/A[k][k]; A[i][k]=g[i][k];}s=0;for(m=0; m<=k-1; m++){s=s+A[k][m]*y[m];}y[k]=(b[k]-s)/A[k][k];}x[n-1]=y[n-1]/A[n-1][n-1];for(k=n-2; k>=0; k--){float sum=0;for(m=k+1; m<n; m++){sum=sum+A[m][k]*x[m];}x[k]=(y[k]-sum)/A[k][k];}for(sq=0; sq<n; sq++){printf("%f ",x[sq]);}return 0;}(3)三对角线性方程组追赶法#include <stdio.h>#include <math.h>#define n 10int main(void){float a[10]={4,4,4,4,4,4,4,4,4,4};float c[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1};float d[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1}; float b[10]={7,5,-13,2,6,-12,14,-4,5,-5}; float x[10]={0};float y[10]={0};float arf[10]={0};float bt[9]={0};arf[0]=a[0];int i;for(i=0;i<n-1;i++){bt[i]=c[i]/arf[i];arf[i+1]=a[i+1]-d[i+1]*bt[i];//printf("%f %f \n",bt[i],arf[i+1]); }y[0]=b[0]/arf[0];//printf("%f\n",y[0]);for(i=1;i<n;i++){y[i]=(b[i]-d[i]*y[i-1])/arf[i];}//printf("%f\n",y[1]);x[n-1]=y[n-1];for(i=n-2;i>=0;i--){x[i]=y[i]-bt[i]*x[i+1]; }for(i=0;i<n;i++)printf("%lf ",x[i]);return 0;}数值结果：(1) 普通的线性方程组①顺序消去法②列主元消去法(2) 对称正定线性方程组平方根法：(3)三对角线性方程组追赶法：对实验计算结果的讨论和分析：(1) 普通的线性方程组①顺序消去法x1~x10的绝对误差：0.000001，-0.000001,0.000001,0,0.000001,0,0.000002,0,0,0x1~x10的相对误差：0.000001，0.000001,-1,0,0.0000005,0,0.00000067,0,0,0误差很小，基本可以忽略。

数值分析(05)高斯消元法

n=length(b); X=zeros(n,1); A的第i行、第i+1到n列元素 X(n)=b(n)/A(n,n); 构成的行向量 for i=n-1:-1:1 X(i)=(b(i)-A(i,i+1:n)* X(i+1:n))/A(i,i); end xn bn / ann
xi (bi

4 7 13 13
解: x4 1 x1 x2 x2 x3 0
3 x3
3 0
x3 2
x1 x2 x2

3
2
求解上三角方程组 Ax=b
for i= n : – 1 : 2
b ( i ) = b ( i ) / A ( i , i );
r ( A) r ( A)方程组Ax b无解（即不相容）。常见是m n，称为超定方程组（又称矛盾方程组）此时，向量b不在A的列空间R( A)之中，原方程组无解，但可求出最小二乘意义下的解 x。即求 x使 || b Ax ||2 2 min
MATLAB实现:
x=A\b
a1n xn b1 a2 n xn b2 ann xn bn
将原方程组 Ax b 化为同解的上三角方组程 Ux g 初等变换 Ax b 同解用增广矩阵表示为
一、三角形方程组的解法
x2 x2 x3 3 x3
x4 5 x4 13 x4 13 x4

4 7 13 13
为求解上三角方程组，从最后一个方程入手，先 a11 x1 a12 x2 .............................. a1n xn b1 解出 xn=bn/ann , 然后按方程由后向前的顺序，从方程 a x2 ............................. a2 n xn b2 22 中依次解出 xn-1,xn-2,…,x ....................................... .. 1。这样就完成了上三角方程组的求解过程。这个过程被称为回代过程其计算步骤如 an1n1 xn1 an1n xn bn 1 下： ann xn bn

数值分析实验报告---高斯消去法 LU分解法

数值分析实验报告---高斯消去法 LU分解法实验一：高斯消去法一、实验目的1. 掌握高斯消去法的原理2. 用高斯消去法解线性方程组3. 分析误差二、实验原理高斯消去法(又称为高斯-约旦消去法)是一种利用矩阵消元的方法，将线性方程组化为改进的阶梯形式，从而解出线性方程组的解的方法。

具体而言，高斯消去法将线性方程组的系数矩阵化为一个上三角矩阵，再利用回带法求解线性方程组的解。

三、实验内容1.1、用高斯消去法解线性方程组在具体实验中，我们将使用高斯消去法来解决下述的线性方程组。

5x+2y+z=102x+6y+2z=14x-y+10z=25为了使用高斯消去法来解这个方程组，首先需要将系数矩阵A进行变换，消除A矩阵中第一列中的下角元素，如下所示:1, 2/5, 1/50, 28/5, 18/50, 0, 49/28接着使用回代法来计算该方程组的解。

回代法的过程是从下往上进行的，具体步骤如下：第三个方程的解：z=49/28；第二个方程的解： y=（14-2z-2x）/6；第一个方程的解： x=（10-2y-z）/5。

1.2、分析误差在使用高斯消去法求解线性方程组时，一般会出现截断误差，导致得到的解与真实解之间存在一些误差。

截断误差的大小和矩阵的维数有关。

为了估计截断误差，我们使用矩阵B来生成误差，在具体实验中，我们将使用下面的矩阵：我们来计算该矩阵的行列式，如果方程组有唯一解，则行列性不为0。

本例中，行列式的值是 -1，因此方程组有唯一解。

然后我们计算真实解和高斯消去法得到的解之间的误差，具体公式如下所示：误差 = 真实解的范数 - 高斯消去法得到的解的范数其中，范数的定义如下：||x||1=max{|xi|}; ||x||2=sqrt{(|x1|^2 + |x2|^2 + ... + |xn|^2)}四、实验步骤1、将高斯消去法的每一个步骤翻译成代码，并保存为一个独立的函数。

2、将代码上传至 Python 交互式环境，并使用高斯消去法来解线性方程组。

解线性代数方程组的直接法之Gauss主元消去法及其C++编程代码

解线性代数方程组的直接法—Gauss 主元消去法最近学习数值分析，为了能够边学习边训练编程，笔者决定对于所学算法一一分析记录，也便于自己以后复习时候的查阅，在这里顺便和大家共享一下我的学习经验，共勉今天主要写的是求解线性方程组的直接法之一Gauss 主元消去法的基本算法步骤和用C++编程实现的具体过程，内容比较简单，代码已经经过测试，可以直接使用。

首先就是Gauss 顺序消去法的基本原理：主元消去法只是在顺序消元法进行前采取选主元措施。

对于Ax=b 这样一非齐次线性方程组，要求解x 列向量，利用线性代数方面的知识，我们可以对(A|b)这样一个增广矩阵做初等行变换，化为矩阵A 下三角为零的形式。

这个时候再进行回代操作就可以计算出来x 列向量啦。

例如：对矩阵Ax=b;111213*********,11,21,31,,1,2,3,A n n n n n n n n n n n n a a a a a a a a a a a a a a a a ----⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ ... ... = ... ... ... ,121...n n b b b b b -⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦ 则对应A 的增广矩阵如下：11112131212223221,11,21,31,1,1,2,3,.........(A |b).........n n n n n n n n n n n n nn b a a a a a a a a b a a a a b a a a a b -----⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第一步：以(A|b)的第一行为基，对的第2行到第n 行进行初等行变换；假设110a ≠，将第1行的111i a a -倍加到第i 行(i=2,3,…,n)，通式为 (2)(1)1111(j 1,2,...,n,n 1)i ij ij j a a a a a =-=+ 注意这里的j 是从1到n+1取值的； a ij 的上标（1）,（2）只是用来表示数据a ij 前一时刻和当前时刻的值经过消元，增广矩阵(A|b)变成下面的形式：(1)(1)(1)(1)(1)11112131(2)(2)(2)(2)222232(2)(2)(2)(2)1,21,31,1(2)(2)(2)(2),2,3, 0...(A |b)......0...0...nnn n n n n n n n nn b a a a a b a a a a a a b a a a b ----⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦到这里，第一步就进行完啦，下来的第2步一直到第n-1步都和上面道理一样，那么进行完n-1步以后，我们得到这样一个增广矩阵：(1)(1)(1)(1)(1)111121,11(2)(2)(2)(2)2222,2(n 1)(n 1)(n 1)1,11,1(n 1)(n 1), 0...(A |b)......00 (00)...0n nnnn n n n n n nn b a a a a b a a a a a b a b ----------⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦嘿嘿，到这里基本就可以大功搞成了。

(完整版)数值计算用C编程实现gauss顺序消去法

用C 编程实现gauss 顺序消去法班级： 09医软（一）班姓名：包优学号： 09713001姓名：曹俊武学号： 09713003姓名：查海松学号： 09713004一、实验名称用c 编程实现gauss 顺序消去法二、执行步骤设有线性方程组Ax = b (1.1)其中A=⎪⎪⎭⎪⎪⎬⎫⎪⎪⎩⎪⎪⎨⎧nn n n n n a a a a a a a a a K K K K KK K 212222111211，x=⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎪⎪⎬⎫n x x x M 21,b=⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎪⎪⎬⎫n b b b M 21 为了清晰起见，以n=4为例来说明消去法的过程，将方程组写成如下的形式⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎪⎪⎬⎫=+++=+++=+++=+++)1(44)1(443)1(432)1(421)1(41)1(34)1(343)1(332)1(321)1(31)1(24)1(243)1(232)1(221)1(21)1(14)1(143)1(132)1(121)1(11b x a x a x a x a b x a x a x a x a b x a x a x a x a b x a x a x a x a (1.2) Gauss 消去法的过程是：第1步假定)1(11a ≠0 ，在方程组（1，2）中保留第一个方程，将第一个方程分别乘 21)1(11)1(21l a a = 、 31)1(11)1(31l a a = 、41)1(11)1(41l a a = ，再分别减第二个、第三个、第四个方程便得到方程组（1,2）的等价方程组⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎪⎪⎬⎫=++=++=++=+++)2(44)2(443)2(432)2(42)2(34)2(343)2(332)2(32)2124)2(243)2(232)2(22)1(14)1(143)1(132)1(121)1(11b x a x a x a b x a x a x a b x a x a x a b x a x a x a x a (1.3)其中)1(11)1()2(ji ij ij a l a a -=(j=2,3,4,5,;i=2,3,4,),)1(11)1()2(b l b b i i i -=(i=2,3,4). 第2 步假定0)2(22≠a ，在方程组（1,3）中保留第一及第二个方程分别乘以 32)2(22)2(32l a a = ，42)2(22)2(42l a a =再分别减第三个、第四个方程便可得到方程组（1,3）的等价方程组⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎪⎪⎬⎫=+=+=++=+++)3(44)3(443)3(43)3(34)3(343)3(33)2124)2(243)2(232)2(22)1(14)1(143)1(132)1(121)1(11b x a x a b x a x a b x a x a x a b x a x a x a x a (1.4) 其中)2(22)2()3(ji ij ij a l a a -=(j=3,4;i=3,4))2(22)2()3(b l b b i i i -=(i=3,4)第3步假定0)3(33≠a ，保留方程组（1,4）的前3个方程，将第三个方程乘 43)3(33)3(43l a a = ，再减去第四个方程便可得到方程组（1,4）的等价方程组 ⎪⎪⎩⎪⎪⎨⎧⎪⎪⎭⎪⎪⎬⎫==+=++=+++)4(44)4(44)3(34)3(343)3(33)2124)2(243)2(232)2(22)1(14)1(143)1(132)1(121)1(11b x a b x a x a b x a x a x a b x a x a x a x a (1.5) 其中)3(3443)3(44)4(44a l a a -=，)3(343)3(4)4(4b l b b -=。

数值分析5(高斯消元法)

a31 a11
b1 b1

an 1 a11
a12 ann
an 1 a11

an 1 a11
消元阶段: 所用除法次数(n-1)+ (n-2)+· · · +1= n(n-1)/2, 所用乘法次数(n-1)*n+ (n-2)* (n-1) +· · · +1*2= (n3-n)/3。
回代阶段:
25/44
13:30
含n个变量的n个方程的一般形式是
a11 a21 a n1 a12 a1n a22 a2 n an 2 ann b1 b2 bn
消去阶段: 使用行初等变换使对角线以下的每个元素都为零。例如,消去第一列对角线以下的元素a21 , a31, · · · , an1。
13:30
求解
a11 x1 a12 x2 a1n xn b1 a x a x a x b 21 1 22 2 2n n 2 an1 x1 an 2 x2 ann xn bn
12/44
13:30
线性方程组的矩阵形式
27/44
计算：xn = bn /ann
a11 x1 a12 x 2 a1n x n b1 a 22 x 2 a 2 n x n b2 a nn x n bn
13:30
高斯消元法由不相等的两部分组成:工作量相对多的消元阶段和工作量相对少的回代阶段。对于大的n, 幂次较低的n相比而言可以忽略不计,最高次幂决定了当n趋近于无穷时的极限形态。换而言之, 对大的n, 低阶项对算法的执行时间的估计没有太大的影响, 当仅需近似估计执行时间时可以忽略不计。

数值分析6(选主元高斯消元法)

x2 1 例2 x1 x2 2
则该方程的精确解为
1 x1 110 4 1.0001, x2 2 x1 0.9999
假设3位十进制浮点数系统,则该方程的浮点形式为
0.100 E 3 x1 + 0.100E 1x2 0.100 E1 0.100 E 1 x1 0.100E 1 x2 0.200 E 1 (0.100E 5 0.100 E 1) x2 0.100 E 5 0.200 E 1 0.100 E 3 x1 + 0.100E 1x2 0.100 E1 0 0.100E 5 x2 0.100 E 5
如果令P =Pn1 Pn2 PA L1 L2 P1 , Ls Pn1 Pn2 Ps1 Fs1 Ps1 Pn2 Pn1 Ln1U LU。
L U Ls是Frobenius矩阵,差别只是Ls 第s列对角线以下元素
是Fs1第s列对角线以下元素经过重新排列得到的。
Matlab 命令: [L,U,P] = lu(A) edit lutx bslashtx
20:22

2 1 2 2 3 3 2 1 2 0 5 0 4 6 19 9
14/27
三对角矩阵分解
f u over , with f
f u( x ) 0 x 1 f (0) a, f (1) b
20:22
k 1
计算次序
1 2
3 4
u11 m 21 A m n1 u12 u22 m n , n 1 u1n u2 n unn
5 6
20:22
10/27
例4 求矩阵的Doolittle分解

高斯消元法_实验报告

- - 华中科技大学数值分析实验报告系、年级研究生院2012级****类别硕士2013年5月6日实验6.1实验要求：根据教材实验6.1做出相应改编：分别使用Gauss 消元、列选主元。

全选主元的方法求解线性方程组，分别比拟三种消元方法的结果和算法的区别，并说明主元的选取在Gauss 消元的中的作用。

问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进展的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。

主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

一般来说书本上采用的列选主元的方法对其线性方程组进展求解的，则我们是否可以选择一种行列都选取主元消去的方法来减小相应的误差呢？全主元消元法和列主元消元法一样都是由高斯消元法演变而来。

只不过选取主元的*围有所加大。

全选主元相对于列选主元的更加复杂化了，因为在运算的过程中导致了元的位置发生了变化，这样我们就不得不追踪每个元的位置。

本次实验就几个问题进展了matlab 实验分析，比拟几种计算方法的优劣性。

实验内容：考虑线性方程组编制一个程序：分别能进展Gauss 消去、列选主元Gauss 消去、全选主元Gauss 消去法进展解线性方程组。

对三种算法所得到的结果进展比拟，分析三种计算方法的准确性。

具体内容：〔1〕取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。

取n=10、n=20计算矩阵的条件数。

分别编写利用matlab 编写运算程序，实现Gauss消去、列选主元消去以及全选主元消去的方法。

比拟三种计算方法的运算结果。

在列选主元的过程中分别采用每步消去过程总选取按模最小或按模尽可能小的元素作为主元或每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

数值分析3 Gauss消去法

§１ Gauss 消去法它是以行初等变换为基础，整理而成的计算机算法。

由它改进和变形得到的高斯选主元消去法及三角分解法，仍然是目前计算机上常用的解低阶稠密矩阵的线性方程组的有效方法。

一、引例：解方程组：⎪⎩⎪⎨⎧-=-+-=+-=+-4336425294321321321x x x x x x x x x用Matlab 先画出空间三个平面，并找出三平面的交点clear,clf,x=[-9 –9 9 9 -9];y=[-9 9 9 –9 -9];z=5-4*x+9*y;plot3(x,y,z),hold on,grid on z=3-2*x+4*y;plot3(x,y,z) z=-4+x-y;plot3(x,y,z)A=[4 –9 2;2 –4 6;-1 1 –3];b=[5 3 –4]';x=A\b; plot3(x(1),x(2),x(3),'r*')解：方程组的增广矩阵是[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-----=431136425294b A消元运算：[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--−−−→−⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡----−−→−-----5.110005.055.00529475.25.225.105.055.0052942313125.025.14142r r r r r r b A回代求解：0.15,1.5/103-=-=x [],5.25.0/)15.0(55.0x 2=-⨯-=95.64/)15.025.295(x 1=⨯+⨯+=。

【注】Gauss 消元的基本思想：用初等行变换方法化一般的方程组为上三角方程组，然后求解；Gauss 消去法步骤：先消元计算，再回代求解二、Gauss 消去法设ｎ阶线性代数方程组为Ax=b (1)其中系数矩阵Ａ(n 阶非奇异矩阵)和右端列向量b 分别为A=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡nn n n n n a a a a a a a a a 212222111211 , b=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡n b b b 21 （一）上三角方程组的回代求解：⎪⎪⎩⎪⎪⎨⎧==++=+++n n nn n n n n b x u b x u x u b x u x u x u 2222211212111 ………………（2）假设u ii ≠０（i=1,2,…，n ），则回代解得：⎪⎩⎪⎨⎧-==∑+=ii ni j j ij i innn n u x u b x u b x /)(/11,2,2,1 --=n n i ………………(3) （二）ＧＡＵＳＳ消元过程：(行变换)记方程组（1）为⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++)1()1(2)1(21)1(1)1(2)1(22)1(221)1(21)1(1)1(12)1(121)1(11nn nn n n n n n n b x a xa x ab x a x a x a b x a x a x a (4)对应的增广矩阵为[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)1()1()1(2)1(1)1(2)1(2)1(22)1(21)1(1)1(1)1(12)1(11)1()1(,n nn n n nn b a a a b a a a b a a a bA…………………(5) 第一轮消元：若)1(11a ≠０，利用)1(11a 消下面的)1(1i a (i=2,3,…n), 消元因子：)1(11)1(11/a a l i i = (i=2,3,…n)[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)2()2()2(2)2(2)2(2)2(22)1(1)1(1)1(12)1(11)2()2(00,n nn n nnb a a b a a b a a a bA………………(6) 其中)1(11)1()2(j i ij ij a l a a -= n j i ,,3,2, =)1(11)1()2(b l b b i i i -= n i ,,3,2 = (7)第二轮消元：若0)2(22≠a ，用第i 行减去第二行的)2(22)2(22/a a l i i =(i=3,4,…n)倍，得到[]⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=)3()3()3(3)3(3)3(3)3(33)2(2)2(2)2(23)2(22)1(1)1(1)1(13)1(12)1(11)3()3(00000,n nnn nn nb a a b a a b a a a b a a a a b A……………(8) 其中n j i b l b b a l a a i i ij i ij ij ,,4,3,,,)2(22)2()3()2(22)2()3( =-=-= …………(9) 已进行k-1轮消元后，得[]⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=)()()()()()()2(2)2(2)2(2)2(22)1(1)1(1)1(1)1(12)1(11)()(,k n k nnk nk k k k knk kk nkn k k K b a a b a a b a a a b a a a a b A……………(10) 第k 轮消元：若)(k kk a ≠0，用第i 行元素减去第k 行元素的)()(/k kk k ik ik a a l = (i=k+1,k+2,…,n)倍，得[]⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=++++++++++++++++)1()1()1(111)1(1)1(11)()()(1)()2(2)2(2)2(12)2(2)2(22)1(1)1(1)1(11)1(1)1(12)1(11)1()1(,k n k nnk nk k k k nk k k k k k kn knk kk k kknk kn k kk K b a a b a a b a a a b a a a a b a a a a a b A…(11) 其中⎪⎩⎪⎨⎧-=-=++,)()()1()()()1(k k ik k i k ik kjik k ij k ij b l b b a l a a ,,,1,n k j i += (12)最后，得[]⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=)()()2(2)2(2)2(22)1(1)1(1)1(12)1(11)()(,n n n nnnnn n b a b a a b a a a bA消元过程：[][][][])()()()()2()2()1()1(,,,,n n k k b A b A b A b A−→−−→−−→−−→−−→− （三）算法公式消元公式：对k=1,2,…,n-1若)(k kk a ≠0 )()(/k kk k ik ik a a l =⎪⎩⎪⎨⎧-=-=++,)()()1()()()1(k k ik k i k ik kj ik k ij k ij b l b b a l a a ,,,1,n k j i +=回代公式：)(n nn a ≠0⎪⎩⎪⎨⎧-==∑+=)(1)()()()(/)(/i ii n i j j i ij i i in nn n n n a x a b x a b x1,2,2,1 --=n n i三、GAUSS 消去法得以实现的前提(可行性分析)1、0a )()2(22(1)11≠n nn a a 或 2、,0,,0,,0,02122221112112221121111≠≠=≠≠A a a a aa a a a a D a a a a a kkk k k kk (13)这是因为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡→⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡)()2(2)2(22)1(1)1(12)1(11212222111211k kk n k kk k k k k a a a a a a a a a a a a a a a 行变换 k i a D n k i ii k ,,2,1,0,0,,,2,1)( =≠≠=则有若令【如】对称正定矩阵、严格对角占优矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡632351214 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-103131021310四、Gauss 消去法的运算量消元过程乘法:[])1(31)()1(211-=-⨯+-∑-=n n k n k n n k 除法:)1(2)(11-=-∑-=n nk n n k 回代过程乘除法总运算量:),1(221+=+++n nnGauss 消去法的总运算量为:3)3(31)1(2)1(2)1(313232n n n n n n n n n n ≈-+=++-+- 对比克莱姆法则的总运算量为(n+1)!【注】对比：当n=10时，克莱姆法则(n+1)!=39916800，而Gauss 消去法，运算４３０次五、算法及例子步１，定义数组A,B,x ；步２，输入A,B,Eps ；步３，对k=1到n-1进行消元：⎪⎩⎪⎨⎧++=-⇐>-⇐n k k j i a b a b b Eps a a a a a a kkk ik i i kk kk kj ik ij ij ,,2,1,/,,/ 步4，nn n n a b x /=；对i=n-1到1回代解ij ni j j iji i a x ab x /)(1∑+=-=步５，输出n i x i ,,2,1, =；步６，结束。

消元法实验报告9

西京学院数学软件实验任务书一、实验问题编写两个程序，分别利用高斯消元法和列主元消去法求解方程组二、分析及计算过程1.顺序消元法源程序：#include "stdio.h"#include "math.h"#include "string.h"void main(){double a[100][100];double x[100];double s,r,e;int m,i,j,k,p,flag=1;printf("请输入矩阵的维数m:\n");scanf("%d",&m);printf("请依次将方程的系数和常数输入：\n");for(i=0;i<m;i++){for(j=0;j<m+1;j++)scanf("%lf",&a[i][j]);}printf("输入的增广矩阵为：\n");for(j=0;j<m;j++){for(i=0;i<m+1;i++)printf("%4.5lf ",a[j][i]);printf("\n");}printf("求的其解：\n");for(k=0;k<m-1;k++){p=k;e=a[k][k];for(i=k+1;i<m;i++){r=a[i][k]/a[k][k];for(j=k;j<m+1;j++)a[i][j]=a[i][j]-r*a[k][j];}}if(flag){x[m-1]=a[m-1][m]/a[m-1][m-1];for(i=m-2;i>=0;i--){s=a[i][m];for(j=i+1;j<m;j++)s=s-a[i][j]*x[j];x[i]=s/a[i][i];}for(i=0;i<m;i++)printf("%15.7lf\n",x[i]);}}例题： 12312312321322323x x x x x x x x x +-=⎧⎪-++=⎨⎪-+=⎩运行结果截图：即：x=1122.高斯列主元消去法源程序：#include "stdio.h" #include "math.h" #include "string.h" void main(){double a[100][100]; double x[100];double s,r,e;int m,i,j,k,p,flag=1;printf("请输入矩阵的维数m:\n");scanf("%d",&m);printf("请依次将方程的系数和常数输入：\n");for(i=0;i<m;i++){for(j=0;j<m+1;j++)scanf("%lf",&a[i][j]);}printf("输入的增广矩阵为：\n");for(j=0;j<m;j++){for(i=0;i<m+1;i++)printf("%4.5lf ",a[j][i]);printf("\n");}printf("求的其解：\n");for(k=0;k<m-1;k++){p=k;e=a[k][k];for(i=k+1;i<m;i++){if(fabs(a[i][k])>e){e=fabs(a[i][k]);p=i;}//选取第一列中最大的for(j=1;j<m+1;j++){s=a[k][j];a[k][j]=a[p][j];a[p][j]=s;}//将第一列中最大的，所在行同第一行交换for(i=k+1;i<m;i++){r=a[i][k]/a[k][k];for(j=k;j<m+1;j++)a[i][j]=a[i][j]-r*a[k][j];}}if(flag)//回带过程{x[m-1]=a[m-1][m]/a[m-1][m-1];for(i=m-2;i>=0;i--){s=a[i][m];for(j=i+1;j<m;j++)s=s-a[i][j]*x[j];x[i]=s/a[i][i];}for(i=0;i<m;i++)printf("%15.7lf\n",x[i]);} }}例题： 1231231232423536x x x x x x x x x -+=⎧⎪--+=⎨⎪++=⎩ 运行截图：即： X=1.111111 0.77777782.5555556。

数值分析实验作业,gauss消去法的数值稳定性分析

实验3.1 Gauss 消去法的数值稳定性试验实验目的:观察和理解Gauss 消元过程中出现小主元(即)(k kka 很小)时引起的方程组解的数值不稳定性。

实验内容：求解方程组b Ax =，其中（1）⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡⨯=11212592.1121-130.6-291.51314.59103.015-1A ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=2178.4617.591b ；（2）⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----=2010151526990999999999.23107102A ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=15019000000000.582b .实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的。

（2）用Gauss 列主元消去法求得L 和U 及解向量421,R x x ∈.（3）用不选主元的Gauss 消去法求得L ~和U ~及解向量421~,~R x x ∈.（4）观察小主元并分析其对计算结果的影响.程序如下:计算矩阵条件数及Gauss 列主元消去法:format longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;k2=cond(A1) %k2为矩阵的条件数；for k=1:n-1a=max(abs(A1(k:n,k))); [p,k]=find(A1==a); B=A1(k,:);c=b1(k);A1(k,:)=A1(p,:);b1(k)=b1(p); A1(p,:)=B;b1(p)=c; if A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); else break end endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1运行结果如下： K2=68.43；⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---⨯=-14929.00202.00893.0011755.04724.00011079.2600118L⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=801.0000231.1835.2001314.5902592.11U 1x =[18.9882;3.3378;-34.747;-33.9865] 不选主元的Gauss 消去法程序：clearformat longengA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2]; n=4;for k=1:n-1A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k);A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); endL1=tril(A1,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A1,0) for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j); endb1(n)=b1(n)/L(n,n); for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1程序运行结果如下：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯=10189.010333.3011168.21033.3700110637.170001~151515L⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡--⨯-⨯-⨯-⨯=-5.000816010637.171091.5210043.101314.59103.0~15151815U ]0;0;0005.1;6848.23[~1=x同理可得2A 对应的系数矩阵条件数及Gauss 列主元消去法求解结果： K2=8.994;⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⨯=1333.04.0001104.0-3.0-0015.0000112-L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=36667.300030.26005.155.2010710U ]0.1;;0.1;0.1;10444.0[152-⨯=-x不选主元的Gauss 消去法结果：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯=1400.0109999.0-001104998.2-5.00013.0-0001~1212L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⨯⨯⨯--=-3667.3000107495.5109987.14003.26100.1010710~121212U ]000145.1;99994.0;000.1;1045.1[~52-⨯-=-x实验4.5 三次样条插值函数的收敛性问题提出：多项式插值不一定收敛的，即插值的节点多，效果不一定就好。

消元法实验报告2

西京学院数学软件实验任务书一、实验目的：了解高斯消元法的基本原理，熟悉高斯顺序消元法和高斯列主元消元法的基本算法。

运用Matlab/C/C++/Java/Maple/Mathematica 等其中一种语言实现高斯消元法。

对于输入的任意矩阵能解其线性方程组。

二、实验基本原理和内容：用高斯消元法求解线性方程组的基本思想是设法消去方程组的系数矩阵A 的主对角线下的元素，而将A x b =化为等价的上三角方程组，然后再回代过程便可以获得方程组的解。

这种解线性方程组的方法，常称为高斯消元法。

内容高斯消元法解方程组： 1230.001 2.000 3.000 1.0001.000 3.712 4.623 2.0002.0001.0725.643 3.000x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦（1）高斯顺序消元法：第一步：消元对1,2,,1k n =-若()0k kk a =则停止计算对1,2,,i k k n =++计算因子()()/k k ik ik kk l a a =；对1,2,,j k k n =++计算(1)()()(1)()()k k k ijij ik kj k k k ii ik k a a l a b b l b ++⎧=-⎪⎨=-⎪⎩；第二步：回代对,1,,1i n n =-()()1ni i iij jj i i i iib a x x a=+-=∑（高斯消元的条件）若A 的所有顺序主子式均不为0，则高斯消元无需换行即可进行到底，且得到唯一解。

若消元过程中允许对增广矩阵进行行交换，则方程组A x b =可用消元法求解的充要条件是A 可逆。

（二）高斯列主元消元法第一步：在A x b =即(1)(1)A x b =的系数矩阵(1)A 的第一列元素中选择一个绝对值最大的元素，不妨设(1)1i a 。

对调(1)1j a 和(1)ij a 及(1)1b 和(1)ib (1,2,,,1)j n i n =≤≤ 。