数值分析上机实验——解线性方程组
- 格式: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分解方法在计算准确性方面与追赶法相当,但在处理较大规模的问题时,计算复杂度较高,追赶法更适合。
【关键字】分析数值分析实验报告二求解线性方程组的直接方法姓名:刘学超日期: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(i=k+1;i n;i++)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");for(i=0;i n;i++)printf("Y%d=%f\n",i+1,x[i]);}结果:四实验收获与教师评语此文档是由网络收集并进行重新排版整理.word可编辑版本!。
《数值分析》实验报告实验编号:实验三课题名称:解线性方程组一、算法介绍1、定义四个函数分别为:DieDai(),Newton(),XianWei(),DuiFen()。
在主函数中输入要选用方法所对应的序号后,用switch语句对函数进行调用。
2、迭代法的主要思想为:保留一个变量在等号的左边,其他都移到等号右边。
3、Newton法的主要算法为:把f(x)在x0附近展开成Taylor级数,取其线性部分,选取一点x0,该点所对应的值为f(x0),对于n=1,2,…,Nmax,按Xn+1=Xn-f(Xn)/f’(Xn)求出Xn+1,并计算f(Xn+1),若|Xn+1-Xn|小于容许误差,则停止计算。
4、弦位法:选定初始值x0,x1,并计算f(x0),f(x1);按迭代公式xn+1=xn-f(xn)(xn-xn-1)/(f(xn)-f(xn-1))计算x2,再求f(x2);如果相邻两次迭代值之差在容许的误差范围内则迭代停止,否则用(x2,f(x2)),(x1,f(x1))分别代替(x1,f(x1)),(x0,f(x0)),重复前两个步骤,直至相邻两次迭代值之差在容许的误差范围内。
5、对分区间法:选取方程的根所在的区间(a,b),取其中点c代入方程中得其值为f(c),如果f(c)与f(a)异号说明方程的根在(a,c)区间中,则令b=c,否则令a=c,如果f(c)的绝对值小于0.0001,则停止运算,否则继续计算,直至f(c)的绝对值小于0.0001。
二、程序代码#include <iostream>#include <iomanip>#include <cmath>using namespace std;double f(double x){x=x*x*x-3*x-1;return x;}void DieDai(){cout<<"迭代序列:\n";double x0=2,x1=0;int i=1;while(fabs(x0-x1)>=0.0001){x1=x0;x0=pow((3*x0+1),1.0/3);cout<<"x"<<i<<"="<<setiosflags(ios::fixed)<<setprecision(4)<<x0<<",f(x"<<i<<")="<<f(x0)<<",误差为:"<<x0*x0*x0-3*x0-1<<endl;i++;}cout<<"方程近似解x*="<<x0<<endl;cout<<"共进行"<<i-1<<"次迭代\n" ;}void Newton(){cout<<"迭代序列:\n";double x0=2,x1=0;int i=1;while(fabs(x0-x1)>=0.0001){x1=x0;x0=x0-f(x0)/(3*x0*x0-3);cout<<"x"<<i<<"="<<setiosflags(ios::fixed)<<setprecision(4)<<x0<<",f(x"<<i<<")="<<f(x0)<<",误差为:"<<x0*x0*x0-3*x0-1<<endl;i++;}cout<<"方程近似解x*="<<x0<<endl;cout<<"共进行"<<i-1<<"次迭代\n" ;}void XianWei(){double x0=1,x1=3,x2=2;cout<<"选定曲线y=f(x)上的两个点P0("<<x0<<","<<f(x0)<<")和P1("<<x1<<","<<f(x1)<<")\n";int i=2;while(fabs(f(x2))>=0.0001){x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));cout<<setiosflags(ios::fixed)<<setprecision(4)<<"当前区间(";cout<<x0<<","<<x1<<"),与x轴交点("<<x2<<","<<f(x2)<<"),误差为:";if(f(x2)*f(x0)<0){cout<<x2*x2*x2-3*x2-1<<endl;x1=x2;}else{cout<<x2*x2*x2-3*x2-1<<endl;x0=x2;}i++;}cout<<"方程近似解x*="<<x2<<endl;cout<<"共进行"<<i-2<<"次迭代\n" ;}void DuiFen(){double a=1,b=3,c;cout<<"f(x)=0的根的存在区间("<<a<<","<<b<<")\n";cout<<"端点函数值f(a)="<<f(a)<<",f(b)="<<f(b)<<endl;int i=2;while(fabs(f(c))>=0.0001){c=(a+b)/2;cout<<setiosflags(ios::fixed)<<setprecision(4)<<"当前区间("<<a<<","<<b<<"),区间中点x="<<c;cout<<",f(x)="<<f(c)<<",误差为:";if(f(c)*f(a)<0){cout<<c*c*c-3*c-1<<endl;b=c;}else{cout<<c*c*c-3*c-1<<endl;a=c;}i++;}cout<<"方程近似解x*="<<c<<endl;cout<<"共进行"<<i-2<<"次迭代\n" ;}void Menu(){int n;cin>>n;switch(n){case 1:cout<<"*** 迭代法***\n";DieDai();Menu();break;case 2:cout<<"*** Newton法***\n";Newton();Menu();break;case 3:cout<<"*** 弦位法***\n";XianWei();Menu();break;case 4:cout<<"*** 对分区间法***\n";DuiFen();Menu();break;case 5:return;}}int main (){cout<<" *****问题:求f(x)=x*x*x-3x-1=0在x0=2附近的实根。
数值分析第一次上机实习报告——线性方程组直接解法一、问题描述设 H n = [h ij ] ∈ R n ×n 是 Hilbert 矩阵, 即11ij h i j =+- 对n = 2,3,4,…13,(a) 取11n n x R ⨯⎛⎫ ⎪=∈ ⎪ ⎪⎝⎭,及n n b H x =,用Gauss 消去法和Cholesky 分解方法来求解n n H y b =,看看误差有多大.(b) 计算条件数:2()n cond H(c) 使用某种正则化方法改善(a)中的结果.二、方法描述1. Gauss 消去法Gauss 消去法一般用于系数矩阵稠密且没有任何特殊结构的线性方程组。
设H =[h ij ],y = (y 1,y 2,…,y n )T . 首先对系数矩阵H n 进行LU 分解,对于k=1,2,…n,交替进行计算:1111),,1,,1(),1,2,,k kj kj kr rj r k ik ik ir rk r kk u h l u j k k n l a l u i k k n u -=-=⎧=-=+⎪⎪⎨⎪=-=++⎪⎩∑∑…… 由此可以得到下三角矩阵L=[l ij ]和上三角矩阵U=[u ij ]. 依次求解方程组Ly=b 和Ux=y ,111,1,2,,1(),,1,,1i i i ir r r n i i ir r r i ii y b l y i n x y u x i n n u -==+⎧=-=⎪⎪⎨⎪=-=-⎪⎩∑∑…… 即可确定最终解。
2. Cholesky 分解法对于系数矩阵对称正定的方程组n n H y b =,存在唯一的对角元素为正数的下三角矩阵L ,使得H=LL T 。
因此,首先对矩阵H n 进行Cholesky 分解,即1122111()1()j jj jj jk k j ij ij ik jk k jj l h l l h l l l -=-=⎧=-⎪⎪⎨⎪=-⎪⎩∑∑ 1,i j n =+… L 的元素求出之后,依次求解方程组Ly=b 和L T x=y ,即1111111(),2,3,i i i ik k k ii b y l y b l y i n l -=⎧=⎪⎪⎨⎪=-=⎪⎩∑… 11(),1,2,n n nn n i i ki k k i nn y x l x y l x i n n l =+⎧=⎪⎪⎨⎪=-=--⎪⎩∑…1 由此求得方程组n n H y b =的解。
数值分析实验报告(一)一、实验名称:求解线性方程组二、实验问题:从2到9的每一个n 值,解n 阶方程组Ax=b 。
这里A 和b 如下定义:a ij =(i+j-1)-1,b i =p(n+i-1)-p(i-1)式中p(x)=(x 2/24)(2+x 2(-7+n 2(14+n(12+3n))))解释发生的现象。
三、计算公式(数学模型):Gauss 消去法是将方程组Ax=b ,用逐次消去未知数的方法,化为等价的三角方程组,再回代求解的方法。
设Ax=b ,A n n R ⨯∈,如果),...,2,1(,0)(n k a k kk=≠,则可通过Gauss 消去法将Ax=b 约化为等价的三角形方程组⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡)()2(2)1(121)()2(2)2(22)1(1)1(12)1(11n n n n nn n n b b b x x x a a a a a a且计算公式为:(a)消元计算(k=1,2,…,n)()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+=-=+=-=+==++n k i a m b b n k j i a m a a n k i a a m k k ik k i k i k kj ik k ij k ijk kn k ik ik ,,1,,1,,1)()()1()()()1()()( (b)回代计算:()⎪⎪⎩⎪⎪⎨⎧-=-==∑+=n i j i ii ji i ij i i in nn n n n n ia x ab x a b x 1)()()()()(1,2,,1四、结构程序设计:1、数据流向图2五、计算程序C++程序//改进的平方根法#include<iostream.h>double p(double x,double n){return x*x/24*(2+x*x*(-7+n*n*(14+n*(12+3*n))));}void main(void){cout<<"*************************线性方程组求解的数值实验题*****************************"<<'\n';double A[9][9]={0};double L[9][9]={0};double B[9]={0};double d[9]={0};double x[9]={0};double y[9]={0};double dd=0;int n,i,j,k;for(n=2;n<10;n++){cout<<n<<"阶方程组Ax=b的解为"<<'\n';for(i=0;i<n;i++){for(j=0;j<n;j++){A[i][j]=1.0/(i+j+1); //赋Ax=b中系数A的初值(A是对称阵)}B[i]=p((double)n+i,(double)n)-p((double)i,(double)n); //赋Ax=b中系数b的初值(调用函数p)}d[0]=A[0][0]; //d[]用来记录A[k][k]for(i=0;i<n;i++) L[i][0]=A[i][0]/d[0]; //LDL'分解中L[i][j]的第一列for(i=1;i<n;i++){for(k=0;k<i;k++) dd+=L[i][k]*L[i][k]*d[k];d[i]=A[i][i]-dd; //LDL'分解中的d[i]dd=0;for(j=i;j<n;j++){for(k=0;k<j;k++) dd+=L[i][k]*d[k]*L[j][k]; //LDL'分解中除第一列外的L[i][j]L[i][j]=(A[i][j]-dd)/d[j];dd=0;}}y[0]=B[0]; //对y的回代过程for(i=1;i<n;i++){for(k=0;k<i;k++) dd+=L[i][k]*y[k];y[i]=B[i]-dd;dd=0;}x[n-1]=y[n-1]/d[n-1]; //对x的回代过程for(i=n-2;i>=0;i--){for(k=i+1;k<n;k++) dd+=L[k][i]*x[k];x[i]=y[i]/d[i]-dd;dd=0;}for(i=0;i<n;i++){cout<<'x'<<i+1<<"="<<x[i]<<'\n';}cout<<'\n';}} //在C++数组中是以开始的,所以i,j,k取值相应变到从0开始运行结果如下:六、结果和讨论:对于2到9阶的线性方程组的解来看,Gauss 消去法有其自己的弊端,使有些结果有误差。
1 / 8数值分析实验六:解线性方程组的迭代法2016113 张威震1 病态线性方程组的求解1.1 问题描述理论的分析表明,求解病态的线性方程组是困难的。
实际情况是否如此,会出现怎样的现象呢?实验内容:考虑方程组Hx=b 的求解,其中系数矩阵H 为Hilbert 矩阵,,,1(),,,1,2,,1i j n n i j H h h i j n i j ⨯===+-这是一个著名的病态问题。
通过首先给定解(例如取为各个分量均为1)再计算出右端b 的办法给出确定的问题。
实验要求:(1)选择问题的维数为6,分别用Gauss 消去法、列主元Gauss 消去法、J 迭代法、GS 迭代法和SOR 迭代法求解方程组,其各自的结果如何?将计算结果与问题的解比较,结论如何?(2)逐步增大问题的维数(至少到100),仍然用上述的方法来解它们,计算的结果如何?计算的结果说明了什么?(3)讨论病态问题求解的算法1.2 算法设计首先编写各种求解方法的函数,Gauss 消去法和列主元高斯消去法使用实验5中编写的函数myGauss.m 即可,Jacobi 迭代法函数文件为myJacobi.m ,GS 迭代法函数文件为myGS.m ,SOR 方法的函数文件为mySOR.m 。
1.3 实验结果1.3.1 不同迭代法球求解方程组的结果比较选择H 为6*6方阵,方程组的精确解为x* = (1, 1, 1, 1, 1, 1)T ,然后用矩阵乘法计算得到b ,再使用Gauss 顺序消去法、Gauss 列主元消去法、Jacobi 迭代法、G-S 迭代法和SOR 方法分别计算得到数值解x1、x2、x3、x4,并计算出各数值解与精确解之间的无穷范数。
Matlab 脚本文件为Experiment6_1.m 。
迭代法的初始解x 0 = (0, 0, 0, 0, 0, 0)T ,收敛准则为||x(k+1)-x(k)||∞<eps=1e-6,SOR方法的松弛因子选择为w=1.3,计算结果如表1。
实验报告哈尔滨工程大学教务处制实验四 解线性方程组一.解线性方程组的基本思想 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.此题为证明题,无程序代码。