数值分析资料报告实验作业,gauss消去法地数值稳定性分析资料报告
Gauss列主元消去法实验
<数值计算方法>实验报告1.实验名称实验2 Gauss 列主元消去法2.实验题目用Gauss 列主元消去法求解线性方程组。
0.0011 2.0002 3.0003 1.0001.0001 3.7122 4.6233 2.0002.0001 1.0722 5.6433 3.000x x x x x x x x x ++=⎧⎪-++=⎨⎪-++=⎩3.实验目的加深自己对Gauss 列主元消去法的理解和认识,并且通过做实验或做练习来加强自己Gauss 列主元消去法的掌握,学会并灵活运用Gauss 列主元消去法来求解方程组。
4.基础理论-------Gauss 列主元消去法1.Gauss 列主元消去法的基本思想是:在进行第k (k=1,2,...,n-1)步消元时,从第k 列的kk a 及以下的各元素中选取绝对值最大的元素,然后通过行变换将它交换到主元素kk a 的位置上,再进行消元。
2.Gauss 列主元消去法的优点:当kk a (k=1,2,...,n-1)的绝对值很小时,用Gauss 列主元消去法来求解方程组时,可以避免所的数值结果产生较大误差或失真。
5.实验环境实验系统:Win 7实验平台:VisualC++语言6.实验过程写出算法→编写程序→计算结果Gauss 列元消去法的算法Input:方程组未知量的个数n;增广矩阵()()1,2,...,T ij A a A A An ==,其中i=1,2,…,n; j=1,2,…,n+1Output:方程组的解x1,x2,…,xn,或失败信息。
1. for i ←1ton-1 do;2. temp ←|ii a |;3. p ←I;4. for j ←i+1 to n do5. if ||ji a >temp then6. p ←j;8. end9. end10. if temp=0 then11. |return False;12. end13. if p ≠I then14. p A ⇔i A ;//i,p 两行交换15. end//列选主元16. for j ←i+1 to n do17.*j ji i A m A -ji m ←/ji ii a a ;18. j A ←*j ji i A m A -;//消元19. end7.实验结果原方程组的解为:X1=-0.490396 , x2=-0.051035 ,x3=0.3675208.附录程序清单#include<iostream.h> #include"stdio.h"#include"math.h"void main ( ){ int n=3,i,j,k,p;doubleA[10][10]={{0.001,2.000,3.000,1.000},{-1.000,3.712,4.623,2.000},{-2.0 00,1.072,5.643,3.000}},temp,m,x[100];for(i=0;i<n;i++){ //选主元temp=fabs(A[i][i]); p=i;for(k=i+1;k<n;k++)if(fabs(A[k][i])>temp){temp=fabs(A[k][i]); p=k;}if(temp==0){ printf("\n无法求解:");return;}if(p!=i)for(j=0;j<n+1;j++){ temp=A[i][j];A[i][j]=A[p][j];A[p][j]=temp;}//消元for(k=i+1;k<n;k++){ m=A[k][i]/A[i][i];for(j=i+1;j<=n;j++)A[k][j]=A[k][j]-m*A[i][j];}}//回代for(i=n-1;i>=0;i--){x[i]=A[i][n];for(j=i+1;j<n;j++)x[i]=x[i]-A[i][j]*x[j];x[i]=x[i]/A[i][i];}printf("\nx=\n");for(i=0;i<n;i++)printf("%f \n",x[i]);}。
数值分析实验作业-gauss消去法的数值稳定性分析说课讲解
数值分析实验作业-g a u s s消去法的数值稳定性分析实验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 三次样条插值函数的收敛性问题提出:多项式插值不一定收敛的,即插值的节点多,效果不一定就好。
高斯消去法的数值稳定性实验报告
end
%½«³£ÊýÏî´ÓAÖзÖÀë
X=A(:,length(A(1,:)));
A=A(:,1:length(A(1,:))-1);
%µÃµ½LÓëU¾ØÕó
fori=1:length(A(:,1))
forj=1:length(A(i,:))
if(i==j)
L(i,j)=1;
U(i,j)=A(i,j);
00-2.83542034035461-0.439277018213440
0000.724838760535495
X2=
5.90739256961544
1.88284104193886
-22.2540722833854
14.5809976298923
(3)beselect.m文件与noselect.m文件
%½«³£ÊýÏî¼Ó½øÈ¥
A(:,length(A(1,:))+1)=Y;
%´ÓÉÏÍùÏÂÑ »·
fori=1:length(A(:,1))
%Ô¤¼ÆËãÿÐеÚÒ»¸öÊý×Ö
fork=i:length(A(:,i))
form=1:i-1
A(k,i) = A(k,i)-A(k,m)*A(m,i);
end
end
end
%¼ÆËãXÖµ
fori=length(X):-1:1
forj=i+1:length(U(i,:))
X(i)=X(i)-U(i,j)*X(j);
end
X(i)=X(i)/U(i,i);
end
end
noselect.m文件:
function[ L,U,X ] = noselect( A,Y )
数值分析实验报告(包括高斯消去、二分法、牛顿迭代法)
for k=1:N
x=(a+b)/2;
fx=feval(f,x);fa=feval(f,a);
if abs((b-a)/2)<e || abs(fx)<e
disp('the number of iterations is');k
f=input('please enter a function:f(x)=');
x0=input('please enter the initial value:x0=');
e=input('please enter error:e=');
N=input('please enter the largest number of iterations:N=');
disp('the approximate solution is');x
disp('f(x) is');fx
disp('the number of iterations is');k
return
else
x0=x;
end
end
end
disp('The maximum number of iterations is reached, stop calculation');
开课学院、实验室:实验时间:2014年1月1日
课程
名称
数值分析基础性实验
实验项目
名称
数值计算算法及实现
数值分析-Gauss消去法
数值分析上机报告1. 考虑方程组⎪⎪⎩⎪⎪⎨⎧-=+++=+++=+++=+++2557.03927.02786.04002.01784.04240.00643.03781.01920.03645.01550.01129.04015.03872..02246.04043.02943.03678.01234.04069.04.3214.3214.3214.321x x x x x x x x x x x x x x x x (1) 用Gauss 消去法解所给方程组(用四位小数计算);(2) 用列主元素消去法解所给方程组并且与(1)比较结果。
1. Matlab 程序>> clearA=input('输入系数矩阵A :');b=input('输入b 向量(按行向量):');B=[A b'];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB ,所以此方程组无解.\n')returnendif RA==RBif RA==nfprintf('请注意:因为RA=RB=%d ,所以此方程组有唯一解.\n',n) X=zeros(n,1);for p=1:n-1for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend%把方程组系数矩阵A 化为同解的上三角矩阵b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);end%从xn 至x1逐个求解上三角方程组elsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')returnendenddisp('方程组的解为:');X运行后输入系数矩阵A:[0.4096 0.1234 0.3678 0.2943;0.2246,0.3872,0.4015,0.1129;0.3645 0.1920 0.3781 0.0643;0.1784 0.4002 0.2786 0.3927]输入b向量(按行向量):[0.4043 0.1550 0.4240 -0.2667]请注意:因为RA=RB=4,所以此方程组有唯一解.方程组的解为:X =-0.2603-1.73242.3280-0.44712.Matlab程序>> clearA=input('输入系数矩阵A:');b=input('输入b向量(按行向量):');B=[A b'];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.\n')endif RA==RBif RA==nfprintf('请注意:因为RA=RB=%d,所以此方程组有唯一解.\n',n)X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1[Y,j]=max(abs(B(p:n,p))); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);endend%把方程组系数矩阵A化为同解的上三角矩阵b=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);end%从xn至x1逐个求解上三角方程组elsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')returnendenddisp('方程组的解为:');X运行结果输入系数矩阵A:[0.4096 0.1234 0.3678 0.2943;0.2246,0.3872,0.4015,0.1129;0.3645 0.1920 0.3781 0.0643;0.1784 0.4002 0.2786 0.3927]输入b向量(按行向量):[0.4043 0.1550 0.4240 -0.2667]请注意:因为RA=RB=4,所以此方程组有唯一解.方程组的解为:X =-1.7324 2.3280 -0.4471。
数值分析实验报告---高斯消去法 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消去法)
《数值分析》实验报告实验编号:实验二课题名称:列主元Gauss消去法一、算法介绍1、输入矩阵的阶数n,方程组的增广矩阵A;2、对k=0,1,…,n-2,循环:选取列中绝对值最大的元素,将主元所在的行的元素保存在数组temp[n+1]中。
若主元为零,则系数矩阵奇异,计算停止;否则,顺序进行。
如果绝对值最大的元素就在矩阵的对角线上,则进行普通高斯消元法的第一大步,否则将方程组系数换行之后再进行普通高斯消元法的第一大步;3、然后利用回代法求解线性方程组。
二、程序代码#include<iostream>#include<cmath>#include<iomanip>using namespace std;int main(){int n=0,k=0,i=0,j=0,h=0,g=0,flag=0,i1,j1;double max=0,m=0;cout<<"***利用列主元Gauss消元法求解线性方程组***"<<endl;cout<<"请输入矩阵的阶数:"<<endl;cin>>n;double a[n][n+1];double t[n+1];double x[n];memset(a,0,sizeof(a));memset(x,0,sizeof(x));cout<<"请输入方程组的增广矩阵:"<<endl;for(i=0;i<n;i++){for(j=0;j<n+1;j++){cin>>a[i][j];}}for(k=0;k<n-1;k++){max=0;j1=0;for(i=k;i<n;i++){if(fabs(a[i][k])>max){max=fabs(a[i][k]);i1=i;j1=k;}}if(max==0){cout<<"该系数矩阵为奇异矩阵,计算停止"<<endl;flag=1;break;}else{cout<<"第"<<j1+1<<"列中绝对值最大的元素是"<<a[i1][j1]<<",在线性方程组的第"<<i1+1<<"行"<<endl;if(i1!=k){for(j=0;j<=n;j++){t[j]=a[i1][j];a[i1][j]=a[k][j];a[k][j]=t[j];}}for(i=k+1;i<=n-1;i++){m=a[i][k]/a[k][k];for(j=k;j<=n;j++)a[i][j]=a[i][j]-m*a[k][j];for(g=0;g<n;g++){for(h=0;h<n+1;h++)cout<<setiosflags(ios::fixed)<<setprecision(2)<<a[g][h]<<" ";cout<<endl;}cout<<endl;}}}if(flag==0){x[n-1]=a[n-1][n]/a[n-1][n-1] ;double sum=0;for(k=n-2;k>=0;k--){sum=0;for(i=n-1;i>=k;i--)sum+=a[k][i]*x[i];x[k]=(a[k][n]-sum)/a[k][k];}cout<<"该线性方程组的解为:"<<endl;for(i=0;i<n;i++)cout<<"x"<<i+1<<"="<<setiosflags(ios::fixed)<<setprecision(2)<<x[i]<<endl;}system("pause");return 0;}三、运算结果截屏四、算法分析列主元Gauss消元法避免了普通高斯消元法中出现的问题:遇到某个主元为零或者当主元绝对值很小时,计算将会停止或求出的结果将与其实际结果相差很远。
数值分析A实验报告
的解。A 及 b 分别为由函数 matrix31 及 array31 所生成的方阵与向量。
function[x,det,flag]=gaussauto(A,b)
[n,m]=size(A);nb=length(b);
if n~=m error('The rows and columns of matrix A must be equal!');
结果分析,由于最大的主元与最小的主元的绝对值没有过于接近于 0,都 是“比较好的一个数”,然而,如果选择 0 能够作为模最小的数的话,则计算 过程将会崩溃。由此可见,虽然矩阵 A 的条件数非常大即显得相当“病态”, 然而,假设不使用 0 元素作为主元素的话,不管以最大模或者最小模作为主元
7
素的话,对于结果都是没有影响的,原因在于不管是哪里一种方法,选取的主 元的绝对值与其他元素相较起来都属于“良态”,当然这个是不考虑 0 元素的 存在的情况。
1.5.4 选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数,重复上述实 验,观察记录并分析实验的结果。
8
我们这里直接引用实验 2 即将出现的 Hilbert 矩阵作为系数矩阵,并给定 其有精确解[1,…,1]’。
我们考察六维度的 Hilbert 系数矩阵,有如下表达,即 H 为
设有
1.0000 0.5000 0.3333 0.2500 0.2000 0.1667 0.5000 0.3333 0.2500 0.2000 0.1667 0.1429 0.3333 0.2500 0.2000 0.1667 0.1429 0.1250 0.2500 0.2000 0.1667 0.1429 0.1250 0.1111 0.2000 0.1667 0.1429 0.1250 0.1111 0.1000 0.1667 0.1429 0.1250 0.1111 0.1000 0.0909
高斯消元法_实验报告
- - 华中科技大学数值分析实验报告系、年级研究生院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消去、列选主元消去以及全选主元消去的方法。
比拟三种计算方法的运算结果。
在列选主元的过程中分别采用每步消去过程总选取按模最小或按模尽可能小的元素作为主元或每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。
Gauss消去法实验报告
程序代码
程序代码说明
function x=gauss(A,b)
n=length(A);
a=[A,b];
for k=1:n-1
maxa=max(abs(a(k:n,k)));
if maxa==0
return;
end
for i=k:n
end
end
x=zeros(n,1);
x(n)=b(n)/a(n,n);
for k=n-1:-1:1
s=b(k);
for j=k+1:n
s=s-a(k,j)*x(j);
end
x(k)=s/a(k,k);
end
x
%调用M文件gauss. m
%形成增广矩阵
%计算乘子
% 对k+1~n项进行消元
%增广矩阵第i行减去第k行的乘子倍目的是将该矩阵中的第k列中a(k,k)以下的元素全部消为零
二、实验原理:
由一般线性方程组在使用Gauss消去法求解时,从求解过程中可以清楚地看到,若 ,必须施以行交换的手续,才能使消去过程继续下去。有时既使 ,但其绝对值很小,由于舍入误差的影响,消去过程也会出现不稳定现象。因此,为使这种不稳定现象发生的可能性减至最小,在施行消去过程时每一步都要选主元素,即要寻找行 ,使
这样就完成了第1步消元。
回代过程:
在最后的一方程中解出 ,得:
再将 的值代入倒数第二个方程,解出 ,依次往上反推,即可求出方程组的解:
其通项为
三、实验内容与步骤:
1、实验内容:依照实验原理编写gauss消去法的程序。
2、实验步骤:首先,在电脑上安装matlab,然后,启动matlab,新建一个M文件。
2020年数值分析实验作业,gauss消去法的数值稳定性分析
作者:败转头作品编号44122544:GL568877444633106633215458 时间:2020.12.13实验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); elsebreakendendL1=tril(A1,0);for i=1:nL1(i,i)=1;endL=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.26000118L ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=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:nL1(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 三次样条插值函数的收敛性问题提出:作者:败转头作品编号44122544:GL568877444633106633215458时间:2020.12.13多项式插值不一定收敛的,即插值的节点多,效果不一定就好。
数值分析实验(Gauss消元法)
end
%这里采用增广矩阵行变换的方式求解
c=n+1;
A(:,c)=b;
%%消元过程
for k=1:n-1
A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c);
end
%%回代结果
x=zeros(length(b),1);
x(n)=A(n,c)/A(n,n);
教师签名
年月日
重庆大学
学生实验报告
实验课程名称数值计算
开课实验室DS1421
学院2009年级土木专业班11班
学生姓名周波学号********
开课时间2010至2011学年第一学期
总成绩
教师签名
课程
名称
数值计算
实验项目
名称
Gauss消元法
实验项目类型
验证
演示
综合
设计
其他
指导
教师
何光辉
成绩
√
实验目的、实验原理
通过高斯消元法的实践环节,达到加深对高斯消元法理解和培养学生程序设计能力的目的。
在Matlab命令窗口中输入Gauss得出结果如下:
>> Gauss
x=
2.0000
1.0000
-1.0000
实验中遇到的问题及解决办法
消元过程要求a(ii)≠0,(i=1,2,…,n-1),回代过程则进一步要求a(nn)≠0。
Matlab求解线性方程的几种命令如下(方程组的一般形式可用矩阵和向量表示成: ,但运用下列方法的前提必须保证所求解的方程为恰定方程,即方程组存在唯一的一组解) :
运用求逆思想 : 或 ;
数值分析实验作业,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 三次样条插值函数的收敛性问题提出:多项式插值不一定收敛的,即插值的节点多,效果不一定就好。
数值分析实验报告
实验五 解线性方程组的直接方法实验5.1 (主元的选取与算法的稳定性) 问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。
主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。
实验内容:考虑线性方程组n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求:(1)取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元,结果如何?(2)现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。
(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。
(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
思考题一:(Vadermonde 矩阵)设⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=∑∑∑∑====n i i n n i i ni i n i i n n n n n n nx x x x b x x x x x x x x x x x x A 002010022222121102001111 ,, 其中,n k k x k ,,1,0,1.01 =+=,(1)对n=2,5,8,计算A 的条件数;随n 增大,矩阵性态如何变化?(2)对n=5,解方程组Ax=b ;设A 的最后一个元素有扰动10-4,再求解Ax=b(3)计算(2)扰动相对误差与解的相对偏差,分析它们与条件数的关系。
高斯消去算法实验报告
高斯消去算法实验报告1. 实验背景高斯消去算法,也称为高斯消元法,是一种用于求解线性方程组的常用方法。
通过进行一系列的行变换,将方程组化简为阶梯矩阵,从而得到方程组的解。
本实验旨在使用高斯消去算法,解决给定的线性方程组。
2. 实验过程2.1 算法原理高斯消去算法的基本思想是通过进行行变换,将线性方程组化简为阶梯矩阵。
具体流程如下:1. 对于每一列,从对角线开始,选取主元(即该列中绝对值最大的元素),并将该主元所在的行与对角线所在的行交换位置。
这样可以避免除法中的误差积累。
2. 通过进行行变换,将主对角线以下的元素全部清零。
具体方法是,对于每一行i,通过消去第i+1行到最后一行的第i列元素,从而将下三角矩阵的元素清零。
3. 倒序遍历每一行,通过行变换,将主对角线以上的元素清零。
具体方法是,消去第i-1行到第1行的第i列元素,从而将上三角矩阵的元素清零。
4. 将矩阵化简为阶梯矩阵。
2.2 实验步骤1. 取得待解线性方程组的系数矩阵A和常数向量b。
2. 将矩阵A和向量b合并为增广矩阵Ab。
3. 通过高斯消去算法,将增广矩阵化简为阶梯矩阵。
4. 根据化简后的阶梯矩阵,求解线性方程组。
3. 实验结果以一个3阶线性方程组为例进行实验,方程组如下:2x + 3y + z = 93x + 2y + 4z = 124x + 3y + 6z = 18按照操作步骤,我们将系数矩阵A和常数向量b合并为增广矩阵Ab:markdownA = [[2, 3, 1],[3, 2, 4],[4, 3, 6]]b = [9, 12, 18]Ab = [[2, 3, 1, 9],[3, 2, 4, 12],[4, 3, 6, 18]]然后,通过高斯消去算法,将增广矩阵Ab化简为阶梯矩阵:markdownAb = [[2, 3, 1, 9],[0, 1.5, 2.5, 6],[0, 0, 0, 0]]根据化简后的阶梯矩阵,我们可以得到方程组的解:x = 1y = 2z = 0因此,该线性方程组的解为x=1,y=2,z=0。
高斯消去实验报告分析
一、实验背景高斯消去法(Gaussian Elimination)是线性代数中解决线性方程组的重要方法之一。
它通过初等行变换将系数矩阵化为行阶梯形矩阵,从而求解线性方程组。
本实验旨在通过编写程序实现顺序高斯消去法和列主元高斯消去法,并对比分析两种方法的计算效率和精度。
二、实验目的1. 理解高斯消去法的原理和步骤。
2. 编写程序实现顺序高斯消去法和列主元高斯消去法。
3. 比较两种方法的计算效率和精度。
4. 分析高斯消去法的适用范围和局限性。
三、实验内容1. 顺序高斯消去法- 实现顺序高斯消去法,将系数矩阵化为行阶梯形矩阵。
- 求解线性方程组,得到未知数的解。
2. 列主元高斯消去法- 实现列主元高斯消去法,将系数矩阵化为行阶梯形矩阵。
- 求解线性方程组,得到未知数的解。
3. 比较两种方法- 对比两种方法的计算效率和精度。
- 分析两种方法的优缺点。
四、实验步骤1. 编写程序- 使用C/C++、Python等编程语言编写顺序高斯消去法和列主元高斯消去法程序。
- 程序应包含以下功能:- 输入线性方程组的系数矩阵和常数项。
- 实现顺序高斯消去法和列主元高斯消去法。
- 输出求解结果。
2. 测试程序- 使用一组已知线性方程组测试程序的正确性。
- 验证求解结果与精确解是否一致。
3. 比较分析- 使用不同的线性方程组测试两种方法的计算效率和精度。
- 比较两种方法的运行时间、解的精度和误差。
五、实验结果与分析1. 顺序高斯消去法- 计算效率:顺序高斯消去法的时间复杂度为O(n^3),其中n为方程组中未知数的个数。
- 精度:顺序高斯消去法可能会受到舍入误差的影响,导致求解结果与精确解存在一定的误差。
2. 列主元高斯消去法- 计算效率:列主元高斯消去法的时间复杂度同样为O(n^3),但通过选择主元,可以减少舍入误差,提高求解精度。
- 精度:列主元高斯消去法具有更高的精度,求解结果与精确解更接近。
3. 比较分析- 顺序高斯消去法和列主元高斯消去法在计算效率上基本相同,但列主元高斯消去法具有更高的精度。
数值分析实验报告
天津工业大学数值分析实验报告电子与信息工程学院例题一用Gauss 消去法来解方程组x 1+x 2+x 3=612x 1−3x 2+3x 3=15−18x 1+3x 2−x 3=−15模型:高斯(Gauss )消去法高斯消去法的基本思想是:先逐次消去变量,将方程组化成同解的上三角形方程组,此过程成为消元过程。
然后按方程相反顺序求解上三角方程组,得到原方程组的解,此过程称为回代过程。
这种方法称为Gauss 消去法。
将方程组改成如下形式a 11 1x 1+a 12 1x 2+······+a 1n 1xn =b 11a 21 1 x 1+a 22 1 x 2+······+a 2n 1 xn =b 2 1 : ∶ ∶ ∶ : ∶ ∶ ∶ a n 1 1x 1+a n 2 1 x 2+······+a nn 1xn =b n1 (2-1) 消元过程:第一步:设a 11 1≠0,记l i 1=a i 1/a 11 1(i =2,3,······,n ),方程中第i 个方程减去第一个方程乘以l i 1,完成第一次消元。
将上述方程式可以改写为a 111x 1+a 12 1x 2+······+a 1n 1xn =b 1 1a 22 2 x 2+······+a 2n 2 xn =b 2 2 ∶ ∶ ∶ ∶ ∶ ∶ a n 22 x 2+······+a nn 2 xn =b n 2(2-2) 其中a ij (2)=a ij (1)−l i 1a 1j (1),b i (2)=b i(1)−l i 1b 1 1(i ,j =2,3·········,n )。
计算方法实验报告(1)----Gauss消去法
计算方法实验报告实验名称:实验(一)Gauss 消去法班级:学生姓名:学号:班级序号:课内序号:指导老师:2018-2019学年第2学期一、实验名称:Gauss消去法二、实验学时: 2学时三、实验目的和要求1、掌握高斯消去法基础原理2、掌握高斯消去法解方程组的步骤3、能用程序语言对Gauss消去法进行编程实现四、实验过程代码及结果1、代码:using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace ConsoleApplication_Gauss{class Program{//回带求值的过程static void CalcX(double[,] a, double[] x, int n){for (int i = n - 1; i >= 0; i--){double sum = 0;for (int j = i + 1; j < n; j++){sum += a[i, j] * x[j];}x[i] = (a[i, n] - sum) / a[i, i];}}//消元的过程static void CalcA(double[,] a, int n){for (int k = 0; k < n - 1; k++){for (int i = k + 1; i < n; i++){//double Lik = a[i, k] / a[k, k];// for (int j = k ; j <= n; j++)for (int j = n; j >= k; j--){a[i, j] = a[i, j] - a[i, k] / a[k, k] * a[k, j];}//a[i, k] = 0;}//Output}}//输出未知数x的值static void Output(double[] x, int n){for (int i = 0; i < n; i++){Console.WriteLine("x[{0}]={1}", i, x[i]);}}static void Output(double[,] a, int n){for (int i = 0; i < n; i++){//string s="";for (int j = 0; j <= n; j++){//s += string.Format("{0,-4}", a[i, j]);Console.Write("{0,6}", a[i, j]);}Console.WriteLine();}}//输入函数,表示输入一串值作为方程组的系数static void Input(double[,] a, int n){for (int i = 0; i <= n - 1; i++){string s = Console.ReadLine();string[] ss = s.Split(' ');for (int j = 0; j <= n; j++){a[i, j] = Convert.ToDouble(ss[j]);}}}static void Main(string[] args){Console.WriteLine("请输入矩阵的维数:");int n =Convert.ToInt32( Console.ReadLine());double[,] a = new double[n,n+1];Console.WriteLine("请输入矩阵的各个元;");Input(a, n);Console.WriteLine("------A(i,j)----------");Output(a, n);CalcA(a, n);Console.WriteLine("------消元之后A(i,j)----------");Output(a, n);double[] x = new double[n];CalcX(a, x, n);Output(x, n);Console.ReadLine();}}}2、结果:…。
数值分析实验,用程序实现Gauss消元法
《数值分析》实验报告实验序号:实验一实验名称:一实验目的:用Microsoft visual c++编写一个功能与Windows 计算器一样的计算器。
二实验内容:(1)算法介绍:计算方法分析:Gauss消元法的基本做法就是把方程组转化成为一个如下图的等价的三角方程组,这个过程叫做消元。
得到三角方程组后,就可以逐个求出Xn,Xn-1,…,X1,这个过程叫回代。
程序代码分析:建立两个数组a和b,通过循环语句将n阶增广矩阵输入进去,通过对列的循环对每一列进行消去未知数,通过n小步n大步把矩阵化简成上三角形矩阵,最后通过迭代法解得方程组得解。
(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次。
数值分析
数值分析实验报告课程:数值计算方法专业:信息与计算科学班级:计算0901 姓名:靳涛(20093213)日期:2010年 12 月 07日数值分析实验报告题目: Gauss 列主元消去法摘要:求解线性方程组的方法很多,主要分为直接法和间接法。
本实验运用直接法的Guass 消去法,并采用选主元的方法对方程组进行求解。
实验目的:1. 学习Gauss 消去法的原理。
2. 了解列主元的意义。
3. 确定什么时候系数阵要选主元数学原理:由于一般线性方程在使用Gauss 消去法求解时,从求解的过程中可以看到,若)1(-k kk a =0,则必须进行行交换,才能使消去过程进行下去。
有的时候即使≠-)1(k kka 0,但是其绝对值非常小,由于机器舍入误差的影响,消去过程也会出现不稳定得现象,导致结果不正确。
因此有必要进行列主元技术,以最大可能的消除这种现象。
这一技术要寻找行r ,使得)1()1(max ||->-=k ikki k rka a并将第r 行和第k 行的元素进行交换,以使得当前的)1(-k kk a 的数值比0要大的多。
这种列主元的消去法的主要步骤如下: 1. 消元过程对k =1,2,…,n -1,进行如下步骤。
1) 选主元,记ikki rk a a >=max ||若||rk a 很小,这说明方程的系数矩阵严重病态,给出警告,提示结果可能不对。
2) 交换增广阵A 的r ,k 两行的元素。
kjrj a a ↔ (j=k,…,n +1)3) 计算消元kkkj ik ij ij a a a a a /-= (i=k+1,…,n ; j =k +1,……,n +1)2. 回代过程对k = n , n -1,…,1,进行如下计算)/(11,∑-=+-=nk j kk j kjn k k a x aa x至此,完成了整个方程组的求解。
程序设计:本实验采用Matlab 的M 文件编写。
Gauss 消去法源程序:cleara=input('输入系数阵:>>\n') b=input('输入列阵b :>>\n') n=length(b); A=[a b]x=zeros(n,1); %%%函数主体for k=1:n-1;%%%是否进行主元选取if abs(A(k,k))<yipusilong;%事先给定的认为有必要选主元的小数yzhuyuan=1;else yzhuyuan=0;endif yzhuyuan; %%%%选主元t=A(k,k); for r=k+1:n;if abs(A(r,k))>abs(t) p=r; else p=k; end end %%%交换元素if p~=k;for q=k:n+1; s=A(k,q);A(k,q)=A(p,q); A(p,q)=s; end end end%%%判断系数矩阵是否奇异或病态非常严重if abs(A(k,k))< yipusilongdisp(‘矩阵奇异,解可能不正确’) end%%%%计算消元,得三角阵for r=k+1:n;m=A(r,k)/A(k,k); for q=k:n+1;A(r,q)=A(r,q)-A(k,q)*m; end end end%%%%求解xx(n)=A(n,n+1)/A(n,n); for k=n-1:-1:1; s=0;for r=k+1:n;s=s+A(k,r)*x(r); endt=(A(k,n+1)-s)x(k)=(A(k,n+1)-s)/A(k,k) end例:求解方程⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡10342212357562z y x ε。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验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); elsebreak end endL1=tril(A1,0); for i=1:nL1(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:nL1(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 三次样条插值函数的收敛性问题提出:多项式插值不一定收敛的,即插值的节点多,效果不一定就好。
对三次样条插值函数又如何呢?理论上证明三次样条插值函数的收敛性是比较困难的,也超过了本课程的容。
通过本实验可以验证这一理论结果.实验容:请按一定的规则分别选择等距或者非等距的插值节点,并不断增加插值节点的个数。
考虑实验4.4中的函数或者选择其他感兴趣的函数,可用Matlab的函数“spline”作此函数的三次样条插值函数。
实验要求:(1)随着节点个数的增加,比较被逼近函数和三次样条差值函数的误差变化情况。
分析所得结果并与拉格朗日插值多项式比较。
(2)三次样条插值函数的思想最早产生于工业部门。
作为工业迎合用的例子,考虑如下例子:某汽车制造商根据三次样条差值函数设计车门曲线,其中一段的数据如下:(3)计算实验4.4的样条插值.程序如下:format shortx1=-1:0.5:1;y1=1./(1+25.*x1.*x1);x2=-1:0.25:1;y2=1./(1+25.*x2.*x2);x3=-1:0.1:1;y3=1./(1+25.*x3.*x3);x4=[-1,-0.82,-0.6,-0.53,-0.34,-0.2,0,0.04,0.2,0.25,0.5,0.8,1];y4=1./(1+25.*x4.*x4);xx=-1:0.01:1;yy1=spline(x1,y1,xx);yy2=spline(x2,y2,xx);yy3=spline(x3,y3,xx);yy4=spline(x4,y4,xx);hold onfplot('1./(1+25.*x.*x)',[-1,1],'m')plot(xx,yy1,'g')plot(xx,yy2,'b')plot(xx,yy3,'k')plot(xx,yy4,'r')hold off %比较被逼近函数与三次样条插值函数图像,直观表现不同插值节点处误差的变化xx=-1:0.2:1;y=1./(1+25.*xx.*xx)%函数在相应节点处的真实值;yy1=spline(x1,y1,xx)y1la=lagrange(x1,y1,xx)yy2=spline(x2,y2,xx)y2la=lagrange(x2,y2,xx)yy3=spline(x3,y3,xx)y3la=lagrange(x3,y3,xx)yy4=spline(x4,y4,xx)y4la=lagrange(x4,y4,xx)其中lagrange函数对应的m文件为:function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;end程序运行结果如下:插值结果比较:y 0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 yy1 0.0385 -0.252 -0.0623 0.3687 0.8024 1 0.8024 0.3687 -0.0623 -0.252 0.0385 y1la 0.0385 -0.3793 -0.1101 0.4005 0.8342 1 0.8342 0.4005 -0.1101 -0.3793 0.0385 yy2 0.0385 0.0551 0.1085 0.1765 0.5342 1 0.5342 0.1765 0.1085 0.0551 0.0385 y2la 0.0385 -0.2587 0.3049 0.0726 0.5636 1 0.5636 0.0726 0.3049 -0.2587 0.0385 yy3 0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 y3la 0.0385 0.0588 0.1 0.2 0.5 1 0.5 0.2 0.1 0.0588 0.0385 yy4 0.0385 0.0587 0.1 0.2016 0.5 1 0.5 0.1989 0.0993 0.0588 0.0385 y4la 0.0385 0.4312 0.1 0.2461 0.5 1 0.5 0.264 0.2387 0.0588 0.0385 x=0:10;y=[0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29];dx0=0.8;dx10=0.2;S=csfit(x,y,dx0,dx10)其中csfit函数的m文件为:function S=csfit(X,Y,dx0,dxn)%Clamped Cubic Spline%Input -X is the 1xn abscissa vector% -Y is the 1xn ordinate vector% -dx0=S'(x0) first derivative boundary condition% -dxn=S'(xn) first derivative boundary condition%Output-S: rows of S are the coefficients, in descending order, for the% cubic interpolantsN=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N));C=H(2:N);U=6*diff(D);%Clamped spline endpoint constraintsB(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1)-dx0);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(dxn-D(N));for k=2:N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);for k=N-2:-1:1M(k+1)=(U(k)-C(k)*M(k+2))/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;for k=0:N-1S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;S(k+1,4)=Y(k+1);endend程序运行结果为:S=。