数值分析上机实验——解线性方程组
迭代法解线性方程组数值分析实验报告

迭代法解线性方程组数值分析实验报告一、实验目的本次实验旨在深入研究和掌握迭代法求解线性方程组的基本原理和方法,并通过数值实验分析其性能和特点。
具体目标包括: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```实验分析:通过高斯消元法,我们成功求解了给定的线性方程组。
数值分析实验报告--解线性方程组的迭代法及其并行算法

disp('请注意:高斯-塞德尔迭代的结果没有达 到给定的精度,并且迭代次数已经超过最大迭 代次数max1,方程组的精确解jX和迭代向量X 如下: ') X=X';jX=jX' end end X=X';D,U,L,jX=jX'
高斯-塞德尔的输入为:
A=[10 2 3;2 10 1;3 1 10]; b=[1;1;2]; X0=[0 0 0]'; X=gsdddy(A,b,X0,inf, 0.001,100) A=[10 2 3;2 10 1;3 1 10]; 请注意:因为对角矩阵 D 非奇异,所以此方程组有解.
0.0301 0.0758 0.1834
8.心得体会:
这已经是第三次实验了, 或多或少我已经对 MATLAB 有了更多的了 解与深入的学习。通过这次实验我了解了雅可比迭代法和高斯- 塞德尔迭代法的基本思想,虽然我们不能熟练编出程序,但还是 能看明白的。运行起来也比较容易,让我跟好的了解迭代法的多 样性,使平常手算的题能得到很好的验证。通过这次实验让我对 MATLAB 又有了更深一层的认识,使我对这门课兴趣也更加浓厚。
运行雅可比迭代程序输入: A=[10
b=[1;1;2];X0=[0 0 0]'; X=jacdd(A,b,X0,inf,0.001,100)
2 3;2 10 1;3 1 10];
结果为:
k= 1 X=
0.1000 k= 2 X= 0.0200 k= 3 X= 0.0400 k= 4 X= 0.0276 k= 5 X= 0.0314 k= 6 X= 0.0294 k= 7 X= 0.0301 k= 8 X= 0.0297
6、 设计思想:先化简,把对角线的项提到左边,其它项
数值分析上机题目

数值分析上机题目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 。
数值分析_线性方程组迭代解法Hilbert矩阵

0.999427 1.000368 1.000068 0.999538 1.003545 0.999834 1.000724 1.002775 0.99582 x 9 0.99888 0.9983 0.99933 0.995748 0.999679
0.996358 0.999739
1.002141 1.001106 1.001204 1.001784 1.005742 1.00315 1.00155 1.00181 1.002206 1.00062
1.000066 0.999977 1.000001 1.000041 0.998534 1.000054 0.999712 0.999029 1.005874 1.000672 1.00175 1.004447
0.996128 0.998886 0.998141 0.994858 10 x 0.993805 0.998912 0.998328 0.998514 0.999227 1.000283 1.000399 1.001125 1.005241 1.001468 1.001894 1.002573 1.006835 1.001551 1.001869 1.002117
数值分析课程上机报告
数值分析第二次上机实习报告
——线性方程组迭代解法
一、问题描述
设 Hn = [hij ] ∈ Rn×n 是 Hilbert 矩阵, 即 hij = 对 n = 2,3,4,…15, 1 i + j −1
1 x ∈ R n×n ,及 bn = H n x ,用 SOR 迭代法和共轭梯度法来求解,并与直 取= 1
四、计算结果及其分析
1. 超松弛迭代法分析 令初始向量 x0=u(1,1,……)T,给定不同的初始向量与松弛因子,计算不同情 况下解的情况。计算结果如下。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验报告哈尔滨工程大学教务处制实验四 解线性方程组一.解线性方程组的基本思想 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.此题为证明题,无程序代码。