数值分析上机实验——解线性方程组
- 格式:doc
- 大小:324.00 KB
- 文档页数:17
迭代法解线性方程组数值分析实验报告一、实验目的本次实验旨在深入研究和掌握迭代法求解线性方程组的基本原理和方法,并通过数值实验分析其性能和特点。
具体目标包括:1、理解迭代法的基本思想和迭代公式的推导过程。
2、掌握雅克比(Jacobi)迭代法、高斯赛德尔(GaussSeidel)迭代法和超松弛(SOR)迭代法的算法实现。
3、通过实验比较不同迭代法在求解不同类型线性方程组时的收敛速度和精度。
4、分析迭代法的收敛性条件和影响收敛速度的因素。
二、实验原理1、线性方程组的一般形式对于线性方程组$Ax = b$,其中$A$ 是$n×n$ 的系数矩阵,$x$ 是$n$ 维未知向量,$b$ 是$n$ 维常向量。
2、迭代法的基本思想迭代法是从一个初始向量$x^{(0)}$出发,按照某种迭代公式逐步生成近似解序列$\{x^{(k)}\}$,当迭代次数$k$ 足够大时,$x^{(k)}$逼近方程组的精确解。
3、雅克比迭代法将系数矩阵$A$ 分解为$A = D L U$,其中$D$ 是对角矩阵,$L$ 和$U$ 分别是下三角矩阵和上三角矩阵。
雅克比迭代公式为:$x^{(k+1)}= D^{-1}(b +(L + U)x^{(k)})$。
4、高斯赛德尔迭代法在雅克比迭代法的基础上,每次计算新的分量时立即使用刚得到的最新值,迭代公式为:$x_i^{(k+1)}=(b_i \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}\sum_{j=i+1}^{n}a_{ij}x_j^{(k)})/a_{ii}$。
5、超松弛迭代法在高斯赛德尔迭代法的基础上引入松弛因子$\omega$,迭代公式为:$x_i^{(k+1)}= x_i^{(k)}+\omega((b_i \sum_{j=1}^{i-1}a_{ij}x_j^{(k+1)}\sum_{j=i}^{n}a_{ij}x_j^{(k)})/ a_{ii} x_i^{(k)})$。
解线性方程组的迭代法数值计算上机实习报告一.综述:考虑用迭代法求解线性方程组,取真解为,初始向量取为零,以范数为度量工具,取误差指标为.其中。
分别完成下面各小题:第六题:编制程序实现Jacobi迭代方法和Gauss-Seidel 方法。
对应不同的停机标准(例如残量,相邻误差,后验误差停机标准),比较迭代次数以及算法停止时的真实误差。
其中残量准则:、相邻误差准则:后验误差停机准则:解:为了结果的可靠性,这里我分别对矩阵阶数为400、2500、10000进行试验,得到对应不同的方法、取不同的停机标准,迭代次数和真实误差的数据如下:分析上面数据可知,对应不同的停机标准,GS方法的迭代次数都近似为J方法的一半,这与理论分析一致。
而且从迭代次数可以看出,在这个例子中,作为停机标准,最强的依次为后验误差,再到残量,再到相邻误差。
第七题:编写程序实现SOR 迭代方法。
以真实误差作为停机标准,数值观测SOR 迭代方法中松弛因子对迭代次数的影响,找到最佳迭代因子的取值。
解:本题中考虑n=50,即对2500阶的矩阵A。
由于我们已经知道要使SOR方法收敛,松弛因子需要位于。
下面来求SOR方法中对应的最佳松弛因子。
应用筛选法的思想,第一次我们取松弛因子,间距为0.05,得到的对应的图像如下,从图中可以看出迭代次数随着的增大,先减小后变大,这与理论相符。
同时可以看出最佳松弛因子.第二次将区间细分为10份,即取,可得下面第二幅图像,从图像中可以看出最佳松弛因子第八题:对于J 方法,GS方法和(带有最佳松弛因子的)SOR 方法,分别绘制误差下降曲线以及残量的下降曲线(采用对数坐标系),绘制(按真实误差)迭代次数与矩阵阶数倒数的关系;解:对于J方法,考虑n=50时,采用相邻误差为迭代的终止条件,误差下降曲线及残量的下降曲线如下:对于GS方法,考虑n=50的时候,采用相邻误差作为迭代的终止条件,所得到的残量和误差的下降曲线如下图:从中可以看出,当相邻误差满足误差指标时,真实误差却并不小于误差指标,而为2.6281e-04。
数值分析实验报告二求解线性方程组的直接方法姓名:刘学超日期:3/28一实验目的1.掌握求解线性方程组的高斯消元法及列主元素法;2.掌握求解线性方程组的克劳特法;3.掌握求解线性方程组的平方根法。
二实验内容1.用高斯消元法求解方程组(精度要求为):2.用克劳特法求解上述方程组(精度要求为)。
3.用平方根法求解上述方程组(精度要求为)。
4.用列主元素法求解方程组(精度要求为):三实验步骤(算法)与结果1用高斯消元法求解方程组(精度要求为):#include stdio.h#define n3 void gauss(double a[n][n],double b[n]){double sum1=0,sum2=0,sum3=0,sum4=0;double l[n][n],z[n],x[n],u[n][n];int i,j,k;for(i=0;i n;i++)l[i][i]=1;for(i=0;i n;i++){for(j=0;j n;j++){if(i=j){for(k=0;k=i-2;k++)sum1+=l[i][k]*u[k][j];u[i][j]=a[i][j]-sum1;}if(i j){for(k=0;k=j-2;k++)sum2+=l[i][k]*u[k][j];l[i][j]=(a[i][j]-sum2)/u[j][j];}}for(k=0;k=i-2;k++)sum3+=l[i][k]*z[k];z[i]=b[i]-sum3;for(i=n-1;i=0;i--){for(k=i;k=n-1;k++)sum4+=u[i][k]*x[k];x[i]=(z[i]-sum4)/u[i][i];}}for(i=0;i n;i++)printf("%.6f",x[i]);}main(){double v[3][3]={{3,-1,2},{-1,2,2},{2,-2,4}};double c[3]={7,-1,0};gauss(v,c);}2用克劳特法求解上述方程组(精度要求为)#include stdio.h#include stdlib.h#include conio.h#define n3 int main(){float u[n][n],l[n][n],d[n]={7,-1,0},x[n];float a[3][3]={{3,-1,2},{-1,2,2},{2,-2,4}};int i,j,k;printf("equations:\n");for(i=0;i n;i++){for(j=0;j n-1;j++)printf("(%f)Y%d+",a[i][j],j+1);printf("(%f)Y%d=%f",a[i][n-1],n,d[i]);printf("\n");}printf("\n");for(j=0;j n;j++)for(i=j;i n;i++)l[i][j]=a[i][j];for(i=0;i n;i++)for(j=i+1;j n;j++)u[i][j]=a[i][j];for(j=1;j n;j++)u[0][j]=u[0][j]/l[0][0];for(k=1;k n;k++){for(j=k;j n;j++)for(i=j;i n;i++)l[i][j]-=l[i][k-1]*u[k-1][j];for(i=k;i n;i++)for(j=i+1;j n;j++)u[i][j]-=l[i][k-1]*u[k-1][j];for(i=k;i n;i++)for(j=i+1;j n;j++)u[k][j]=u[k][j]/l[k][k];}d[0]=d[0]/l[0][0];for(k=0;k 2;k++){for(i=k+1;i n;i++)d[i]-=d[k]*l[i][k];d[k+1]/=l[k+1][k+1];}for(i=0;i n;i++)x[i]=d[i];for(k=n-2;k 2-n;k--)for(i=k;i-1;i--)x[i]-=x[k+1]*u[i][k+1];for(j=0;j n;j++)for(i=j;i n;i++)printf("l[%d][%d]=%f\n",i+1,j+1,l[i][j]);printf("\n");for(i=0;i n;i++)for(j=i+1;j n;j++)printf("u[%d][%d]=%f\n",i+1,j+1,u[i][j]);printf("\n");for(i=0;i n;i++)printf("d%d=%f\n",i+1,d[i]);printf("\n");printf("the result is:\n");for(i=0;i n;i++)printf("Y%d=%f\n",i+1,x[i]);getch();}结果:3用平方根法求解上述方程组(精度要求为)#include stdio.h#define n3 void gauss(double a[n][n],double b[n]) {double sum1=0,sum2=0,sum3=0,sum4=0;double l[n][n],z[n],x[n],u[n][n];int i,j,k;for(i=0;i n;i++)l[i][i]=1;for(i=0;i n;i++){for(j=0;j n;j++){if(i==j){for(k=0;k=i-2;k++)sum1+=pow(l[i][k],2);l[i][j]=sqrt(a[i][i]-sum1);}if(i j){for(k=0;k=j-2;k++)sum2+=l[i][k]*u[k][j];l[i][j]=(a[i][j]-sum2)/l[j][j];}}for(k=0;k=i-2;k++)sum3+=l[i][k]*z[k];z[i]=(b[i]-sum3)/l[i][i];for(i=n-1;i=0;i--){for(k=i;k=n-1;k++)sum4+=l[k][i]*x[k];x[i]=(z[i]-sum4)/l[i][i];}}for(i=0;i n;i++)printf("%.6f",x[i]);}main(){double v[3][3]={{3,-1,2},{-1,2,2},{2,-2,4}};double c[3]={7,-1,0};gauss(v,c);}结果:4用列主元素法求解方程组(精度要求为):#include stdio.h#include math.h#define n3 int main(){float u[n][n],l[n][n],d[n]={7,-1,0},x[n];float a[n][n]={3,-1,2,-1,2,-2,2,-2,4};int i,j,k;printf("equations:\n");for(i=0;i n;i++){for(j=0;j n-1;j++)printf("(%f)Y%d+",a[i][j],j+1);printf("(%f)Y%d=%f",a[i][n-1],n,d[i]);printf("\n");}printf("\n");for(i=0;i n;i++)for(j=0;j n;j++)l[i][j]=a[i][j];for(i=0;i n;i++)for(j=0;j n;j++)u[i][j]=a[i][j];l[0][0]=sqrt(l[0][0]);u[0][0]=sqrt(u[0][0]);for(i=1;i n;i++)l[i][0]/=u[0][0];for(j=1;j n;j++)u[0][j]/=l[0][0];for(k=1;k 3;k++){for(j=0;j k;j++)l[k][k]-=pow(l[k][j],2);l[k][k]=sqrt(l[k][k]);for(j=0;j k;j++)l[i][k]-=l[i][j]*l[k][j];for(i=k+1;i n;i++)for(j=0;j k;j++)l[i][k]/=l[k][k];}d[0]=d[0]/l[0][0];for(k=0;k 2;k++){for(i=k+1;i n;i++)d[i]-=d[k]*l[i][k];d[k+1]/=l[k+1][k+1];}for(i=0;i n;i++)for(j=0;j n;j++)u[i][j]=l[j][i];for(k=n-1;k 1-n;k--){x[k]=d[k]/u[k][k];for(i=k-1;i-1;i--)d[i]=d[i]-u[i][k]*x[k];}for(j=0;j n;j++){for(i=j;i n;i++)printf("l[%d][%d]=%f\n",i+1,j+1,l[i][j]);}printf("\n");for(i=0;i n;i++){for(j=i;j n;j++)printf("u[%d][%d]=%f\n",i+1,j+1,u[i][j]);}printf("\n");printf("the result is:\n");printf("Y%d=%f\n",i+1,x[i]);}结果:四实验收获与教师评语。
实验题目:线性代数求解方程组一、实验目的1. 理解线性代数中方程组的求解方法。
2. 掌握利用计算机求解线性方程组的算法。
3. 熟悉数学软件(如MATLAB、Python等)在数学问题中的应用。
二、实验内容本次实验主要利用数学软件求解线性方程组。
线性方程组是线性代数中的一个基本问题,其求解方法有很多种,如高斯消元法、矩阵求逆法等。
本实验以高斯消元法为例,利用MATLAB软件求解线性方程组。
三、实验步骤1. 编写高斯消元法算法程序。
2. 输入方程组的系数矩阵和常数项。
3. 调用程序求解方程组。
4. 输出解向量。
四、实验代码及分析1. 高斯消元法算法程序```matlabfunction x = gaussElimination(A, b)[n, m] = size(A);assert(n == m, 'The matrix A must be square.');assert(n == length(b), 'The length of b must be equal to the number of rows in A.');% 初始化解向量x = zeros(n, 1);% 高斯消元for i = 1:n-1% 寻找最大元素[~, maxIdx] = max(abs(A(i:n, i)));maxIdx = maxIdx + i - 1;% 交换行A([i, maxIdx], :) = A([maxIdx, i], :);b([i, maxIdx]) = b([maxIdx, i]);% 消元for j = i+1:nfactor = A(j, i) / A(i, i);A(j, i:n) = A(j, i:n) - factor A(i, i:n); b(j) = b(j) - factor b(i);endend% 回代求解for i = n:-1:1x(i) = (b(i) - A(i, i+1:n) x(i+1:n)) / A(i, i); endend```2. 输入方程组的系数矩阵和常数项```matlabA = [2, 1, -1; 1, 2, 1; -1, 1, 2];b = [8; 5; 2];```3. 调用程序求解方程组```matlabx = gaussElimination(A, b);```4. 输出解向量```matlabdisp('解向量为:');disp(x);```五、实验结果与分析实验结果:```解向量为:2-13```实验分析:通过高斯消元法,我们成功求解了给定的线性方程组。
数值分析上机题目4(总21页) --本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--实验一实验项目:共轭梯度法求解对称正定的线性方程组 实验内容:用共轭梯度法求解下面方程组(1) 123421003131020141100155x x x x -⎛⎫⎛⎫⎛⎫ ⎪ ⎪ ⎪--- ⎪ ⎪ ⎪=⎪ ⎪ ⎪-- ⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ 迭代20次或满足()(1)1110k k x x --∞-<时停止计算。
编制程序:储存m 文件function [x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r; k=0;while rho>10^(-11) & k<1000 k=k+1; if k==1 p=r; elsebeta=rho/rho1; p=r+beta*p; end w=A*p;alpha=rho/(p'*w); x=x+alpha*p; r=r-alpha*w; rho1=rho;rho=r'*r; end运行程序: clear,clcA=[2 -1 0 0;-1 3 -1 0;0 -1 4 -1;0 0 -1 5]; b=[3 -2 1 5]'; [x,k]=CGmethod(A,b)运行结果: x =(2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵, A[i,i]=4,A[i,i-1]=A[i-1,i]=-1,i=2,3,..,n b[1]=3, b[n]=3, b[i]=2,i=2,3,…,n-1迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。
编制程序:储存m 文件function [x,k]=CGmethod_1(A,b) n=length(A);x(1:n,1)=0;r=b-A*x;r1=r; k=0;while norm(r1,1)>10^(-7)&k<10^4 k=k+1; if k==1 p=r; elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p; end r=r1; w=A*p;alpha=(r'*r)/(p'*w); x=x+alpha*p; r1=r-alpha*w; end运行程序: clear,clc n=1000; A=hilb(n); b=sum(A')';[x,k]=CGmethod(A,b)实验二1、 实验目的:用复化Simpson 方法、自适应复化梯形方法和Romberg 方法求数值积分。
(完整word版)迭代法解线性方程组-数值分析实验报告编辑整理:尊敬的读者朋友们:这里是精品文档编辑中心,本文档内容是由我和我的同事精心编辑整理后发布的,发布之前我们对文中内容进行仔细校对,但是难免会有疏漏的地方,但是任然希望((完整word版)迭代法解线性方程组-数值分析实验报告)的内容能够给您的工作和学习带来便利。
同时也真诚的希望收到您的建议和反馈,这将是我们进步的源泉,前进的动力。
本文可编辑可修改,如果觉得对您有帮助请收藏以便随时查阅,最后祝您生活愉快业绩进步,以下为(完整word版)迭代法解线性方程组-数值分析实验报告的全部内容。
数学与计算科学学院《数值分析》课程设计题目:迭代法解线性方程组专业:信息与计算科学学号: 1309302—24姓名:谭孜指导教师:郭兵成绩:二零一六年六月二十日一、前言:(目的和意义)1.实验目的①掌握用迭代法求解线性方程组的基本思想和步骤.②了解雅可比迭代法,高斯—赛德尔法和松弛法在求解方程组过程中的优缺点。
2。
实验意义迭代法是用某种极限过程去逐步逼近线性方程组精确解的方法,它是解高阶稀疏方程组的重要方法。
迭代法的基本思想是用逐次逼近的方法求解线性方程组。
比较雅可比迭代法,高斯—赛德尔迭代方法和松弛法,举例子说明每种方法的试用范围和优缺点并进行比较.二、数学原理:设有方程组b Ax = …① 将其转化为等价的,便于迭代的形式f Bx x += …② (这种转化总能实现,如令b f A I B =-=,), 并由此构造迭代公式f Bx x k k +=+)()1( …③ 式中B 称为迭代矩阵,f 称为迭代向量。
对任意的初始向量)0(x ,由式③可求得向量序列∞0)(}{k x ,若*)(lim x x k k =∞→,则*x 就是方程①或方程②的解。
此时迭代公式②是收敛的,否则称为发散的。
构造的迭代公式③是否收敛,取决于迭代矩阵B 的性 1。
雅可比迭代法基本原理设有方程组),,3,2,1(1n i b x a j j nj ij ==∑= …①矩阵形式为b Ax =,设系数矩阵A 为非奇异矩阵,且),,3,2,1(,0n i a ii =≠从式①中第i 个方程中解出x,得其等价形式)(111j nj j ij ii i x a b a x ∑≠=-= …②取初始向量),,,()0()0(2)0(1)0(n x x x x =,对式②应用迭代法,可建立相应的迭代公式: )(111)()1(∑≠=++-=nj j i k j ij ii k ib x a a x…③ 也可记为矩阵形式:J x J k F B x k +==)()1( …④ 若将系数矩阵A 分解为A=D —L-U ,⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=--=--00000000000000111211212211212222111211n n n nn n nn nn n n n n a a a a a a a a a a a a a a a a a a U L D A式中⎪⎪⎪⎪⎪⎭⎫⎝⎛=nn a a a D2211,⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=-0000121323121nn n n a a a a a a L ,⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=-0000122311312n n n n a a a a a a U 。
数值分析实验报告-清华大学--线性代数方程组的数值解法(总15页)--本页仅作为文档封面,使用时请直接删除即可----内页可以根据需求调整合适字体及大小--线性代数方程组的数值解法实验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)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A);Ab=[A b];r=input('选主元方式(0:自动;1:手动),r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元,请输入第',num2str(i),'列所选元素所处的行数:']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =×103Cond(A,2) = ×103Cond(A,inf) =×103程序自动选择主元(列主元)a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=0b.计算结果x=[1,1,1,1,1,1,1,1,1,1]T(2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2…(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[, , , , , , , , , ]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k 步选择主元处于第k+1行) 最终计算得x=[1,1,1,1,1,1,1,1,1,1]T(3)n=20,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2…(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[,,,,,,,,,,,,,,,,,,,]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k步选择主元处于第k+1行)最终计算得x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
数值分析上机实践报告一、实验目的本实验的目的是通过编写数值分析程序,掌握解决数学问题的数值计算方法,并通过实际应用来检验其有效性和准确性。
具体包括以下几个方面的内容:1.掌握二分法和牛顿迭代法的基本原理和实现方法;2.熟悉利用矩阵的LU分解和追赶法解线性方程组的过程;3.通过具体的实例应用,比较不同方法的计算效果和精度。
二、实验内容本实验分为三个部分,每个部分包括一个具体的数学问题和相应的数值计算方法。
1.问题一:求方程f(x)=x^3-5x^2+10x-80=0的近似解。
在问题一中,我们通过二分法和牛顿迭代法来求解方程的近似解,并比较两种方法的精度和收敛速度。
2.问题二:用LU分解解线性方程组。
问题二中,我们通过矩阵的LU分解方法解线性方程组Ax=b,然后和直接用追赶法解线性方程组进行对比,验证LU分解的有效性和准确性。
三、实验结果及分析1.问题一的结果分析:通过二分法和牛顿迭代法求解方程f(x)=x^3-5x^2+10x-80=0的近似解,得到的结果如下:从结果来看,两种方法得到的近似解均与真实解x≈5非常接近。
但是,通过比较可以发现,牛顿迭代法的计算速度比二分法更快,迭代的次数更少。
因此,在需要高精度近似解的情况下,牛顿迭代法是一个更好的选择。
2.问题二的结果分析:通过LU分解和追赶法解线性方程组Ax=b,得到的结果如下:-用LU分解解线性方程组得到的结果为x1≈1.0,x2≈2.0,x3≈3.0;-用追赶法解线性方程组得到的结果为x1≈1.0,x2≈2.0,x3≈3.0。
从结果来看,两种方法得到的结果完全一致,而且与真实解非常接近。
这表明LU分解方法和追赶法均可以有效地解决线性方程组问题。
但是,在实际应用中,当方程组规模较大时,LU分解方法的计算复杂度较高,因此追赶法更加适用。
四、实验总结通过本实验,我掌握了二分法和牛顿迭代法以及LU分解和追赶法的基本原理和实现方法。
通过具体的数学问题实例应用,我比较了不同方法的计算效果和精度,得出以下结论:1.在求解函数的近似解时,牛顿迭代法相对于二分法具有更快的收敛速度和更高的计算精度;2.在解决线性方程组问题时,LU分解方法在计算准确性方面与追赶法相当,但在处理较大规模的问题时,计算复杂度较高,追赶法更适合。
实验报告哈尔滨工程大学教务处制实验四 解线性方程组一.解线性方程组的基本思想 1.直接三角分解法:将系数矩阵A 转变成等价两个矩阵L 和U 的乘积 ,其中L 和U 分别是下三角和上三角矩阵。
当A 的所有顺序主子式都不为0时,矩阵A 可以分解为A=LU ,且分解唯一。
其中L 是单位下三角矩阵,U 是上三角矩阵。
2.平方根法:如果矩阵A 为n 阶对称正定矩阵,则存在一个对角元素为正数的下三角实矩阵L ,使得:A=LL^T 。
当限定L 的对角元素为正时,这种分解是唯一的,称为平方根法(Cholesky )分解。
3.追赶法:设系数矩阵为三对角矩阵1122233111000000000000000n n n nn b c a b c a b A a b c a b ---⎛⎫ ⎪ ⎪ ⎪=⎪ ⎪⎪⎪ ⎪⎝⎭则方程组Ax=f 称为三对角方程组。
设矩阵A 非奇异,A 有Crout 分解A=LU ,其中L 为下三角矩阵,U 为单位上三角矩阵,记1122233110000100000001000000100,00000000000001n n nn b L U γαβγββγβ--⎛⎫⎛⎫ ⎪⎪ ⎪ ⎪ ⎪ ⎪∂==⎪⎪ ⎪⎪ ⎪ ⎪⎪ ⎪⎪⎪∂⎝⎭⎝⎭ 可先依次求出L ,U 中的元素后,令Ux=y ,先求解下三角方程组Ly=f 得出y ,再求解上三角方程组Ux=y 。
4.雅克比迭代法:首先将方程组中的系数矩阵A 分解成三部分,即:A = L+D+U ,如图1所示,其中D 为对角阵,L 为下三角矩阵,U 为上三角矩阵。
之后确定迭代格式,X )1(+k = BX )(k +f ,如图2所示,其中B 称为迭代矩阵,雅克比迭代法中一般记为J 。
(k = 0,1,......)再选取初始迭代向量X )0(,开始逐次迭代。
5.超松弛迭代法(SOR )它是在GS 法基础上为提高收敛速度,采用加权平均而得到的新算法。
选取分裂矩阵M 为带参数的下三角矩阵M =ω1(D -L ω), 其中ω>0 为可选择的松弛因子,一般当1<ω<2时称为超松弛。
二.实验题目及实验目的1.(第五章习题8)用直接三角分解(杜利特尔(Doolittle )分解)求线性方程组141x +251x +361x = 9, 131x +241x +351x = 8,121x + 2x +32x = 8 的解。
2.(第五章习题9)用追赶法解三对角方程组Ax=b ,其中A=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--------2100012100012100012100012,b=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛00001. 3.(第五章习题10)用改进的平方根法解线性方程组⎪⎪⎪⎭⎫ ⎝⎛---131321112⎪⎪⎪⎭⎫ ⎝⎛321x x x = ⎪⎪⎪⎭⎫ ⎝⎛654 4.(第六章习题7)用SOR 方法解线性方程组(分别取松弛因子ω=1.03,ω=1,ω=1.1)41x - 2x = 1, -1x +42x - 3x = 4,-2x +43x = -3.精确解x *=(21,1,-21)T .要求当)(*k x x -∞<5×106-时迭代终止,并且对每一个ω值确定迭代次数.5.(第六章习题8)用SOR 方法解线性方程组(取ω=0.9) 51x -22x + 3x = -12, -1x +42x - 23x = 20, 21x -32x +103x = 3.要求当)()1(k k x x -+∞<104-时迭代终止.6.(第六章习题9)设有线性方程组Ax=b ,其中A 为对称正定阵,迭代公式)()1(k k x x =++ω(b- A )(k x ),k=0,1,2…,试证明当0<ω<β2时上述迭代法收敛(其中0<α≤λ(A)≤β). 7.(第六章计算实习题1)给出线性方程组H n x=b ,其中系数矩阵H n 为希尔伯特矩阵:H n x=(h ij )∈R n n ⨯, h ij =11-+j i ,i ,j=1,2,…,n.假设x *=(1,1,…,1)T ∈R n ,b= H n x *.若取n=6,8,10,分别雅克比迭代法及SOR 迭代(ω=1,1.25,1.5)求解.比较计算结果. 三.实验手段:指操作环境和平台:win7系统下MATLAB R2009a程序语言:一种类似C 语言的程序语言,但比C 语言要宽松得多,非常方便。
四.程序1.①直接三角分解(文件ZJsanjiao.m ) function x=ZJsanjiao(A,b) [m,n]=size(A); [l u]=lu(A); s=inv(l)*[A,b]; x=ones(m,1);for i=m:-1:1h=s(i,m+1);for j=m:-1:1;if j~=ih=h-x(j)*s(i,j);endendx(i)=h/s(i,i);end②控制台输入代码:>> A=[1/4,1/5,1/6;1/3,1/4,1/5;1/2,1,2]; >> b=[9;8;8];>> x=ZJsanjiao(A,b)2.①追赶法(文件ZG_SDJ.m)function x=ZG_SDJ(a,b,c,f)%aÊǶԽÇÏßÔªËØ%bÊǶԽÇÏßÉÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö%cÊǶԽÇÏßÏ·½µÄÔªËØ£¬¸öÊý±ÈaÉÙÒ»¸ö%fÊdz£ÊýÏîbN=length(a);b=[b,0];c=[0,c];a1=zeros(N,1);b1=zeros(N,1);y=zeros(N,1);x=zeros(N,1);a1(1)=a(1);b1(1)=b(1)/a1(1);y(1)=f(1)/a1(1);for j1=2:Na1(j1)=a(j1)-c(j1)*b1(j1-1);b1(j1)=b(j1)/a1(j1);temp1=f(j1)-c(j1)*y(j1-1);y(j1)=temp1/a1(j1);endj1=N;x(j1)=y(j1);for j1=N-1:-1:1x(j1)=y(j1)-b1(j1)*x(j1+1);end②控制台输入代码:>> a=[2 2 2 2 2];>> b=[-1 -1 -1 -1];>> c=[-1 -1 -1 -1];>> f=[1;0;0;0;0];>> x=ZG_SDJ(a,b,c,f)3.①改进的平方根法(文件GJPFG.m)function GJPFG(A,b)n=length(b);% nΪÁÐά£»% LDL'·Ö½â£»d(1)=A(1,1);for i=2:nfor j=1:i-1sum1=0;for k=1:j-1sum1=sum1+t(i,k)*l(j,k);endt(i,j)=A(i,j)-sum1;l(i,j)=t(i,j)/d(j);endsum2=0;for k=1:i-1sum2=sum2+t(i,k)*l(i,k);endd(i)=A(i,i)-sum2;endfor i=1:nl(i,i)=1;enddisp('µ¥Î»ÏÂÈý½Ç¾ØÕóLΪ£º'); %½â³öµ¥Î»ÏÂÈý½Ç¾ØÕóL£»ldisp('¶Ô½Ç¾ØÕóDΪ£º'); %½â³ö¶Ô½Ç¾ØÕóD£»d%ÓÉLDL'x=bÇó½âx£»%ÓÉLy=b£¬Çóy£»%ÓÉL'x=inv£¨D£©y£¬Çó½âx£»y(1)=b(1);for i=2:nsum3=0;for k=1:i-1sum3=sum3+l(i,k)*y(k);endy(i)=b(i)-sum3;endx(n)=y(n)/d(n);for i=n-1:-1:1sum4=0;for k=i+1:nsum4=sum4+l(k,i)*x(k);endx(i)=(y(i)/d(i))-sum4;enddisp('ÓÉLy=bÇó½âyµÃ£º');ydisp('Ax=bµÄ½âxΪ£º');x②控制台输入代码:>> A=[2 -1 1;-1 -2 3;1 3 1];>> b=[4;5;6];>> GJPFG(A,b)4.①SOR方法(文件SOR_1.m)function SOR_1(A,b,x0,x_a,w)%x_aΪ¾«È·½âif(w<=0 || w>=2)error('²ÎÊý·¶Î§´íÎó');return;endeps=5.0e-6;D=diag(diag(A)); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄÉÏÈý½ÇÕóB=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1; %µü´ú´ÎÊýwhile norm(x_a-x)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=200)disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡');return;endenddisp('Ax=bµÄ½âΪ£º');xdisp('µü´ú´ÎÊýΪ£º');n②控制台输入代码:>> A=[4 -1 0;-1 4 -1;0 -1 4];>> b=[1;4;-3];>> x0=[0;0;0];>> x_a=[0.5;1;-0.5];>> w=1.03;>> SOR_1(A,b,x0,x_a,w)>> w=1;>> SOR_1(A,b,x0,x_a,w)>> w=1.1;>> SOR_1(A,b,x0,x_a,w)5.①SOR方法(文件SOR_2.m)function SOR_2(A,b,x0,w,eps)if(w<=0 || w>=2)error('²ÎÊý·¶Î§´íÎó');return;endD=diag(diag(A)); %ÇóAµÄ¶Ô½Ç¾ØÕóL=-tril(A,-1); %ÇóAµÄÏÂÈý½ÇÕóU=-triu(A,1); %ÇóAµÄÉÏÈý½ÇÕóB=inv(D-L*w)*((1-w)*D+w*U);f=w*inv((D-L*w))*b;x=B*x0+f;n=1; %µü´ú´ÎÊýwhile norm(x-x0)>=epsx0=x;x =B*x0+f;n=n+1;if(n>=200)disp('Warning: µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²£¡');return;endenddisp('Ax=bµÄ½âΪ£º');xdisp('µü´ú´ÎÊýΪ£º');n②控制台输入代码:>> A=[5 2 1;-1 4 2;2 -3 10];>> b=[-12;20;3];>> x0=[0;0;0];>> w=0.9;>> eps=10e-4;>> SOR_2(A,b,x0,w,eps)6.此题为证明题,无程序代码。