雅可比迭代实验报告

第1篇一、实验目的1. 理解雅各比迭代法的原理和应用。

2. 掌握雅各比迭代法的计算步骤和实现方法。

3. 通过实验验证雅各比迭代法在求解线性方程组中的有效性和收敛性。

二、实验原理雅各比迭代法是一种求解线性方程组的迭代方法。

对于形如Ax=b的线性方程组，其中A是n×n的系数矩阵，x是n维未知向量，b是n维常数向量，雅各比迭代法的基本思想是将方程组Ax=b转化为一系列的简单方程进行迭代求解。

设A为对角占优矩阵，则雅各比迭代法的迭代公式为：x_{k+1} = (D - L)^{-1}(b - Ux_k)其中，D是A的对角矩阵，L是A的非对角元素中下三角矩阵，U是A的非对角元素中上三角矩阵。

三、实验内容1. 准备实验环境：安装MATLAB软件，创建实验文件夹。

2. 编写实验程序：（1）定义系数矩阵A和常数向量b。

（2）计算对角矩阵D、下三角矩阵L和上三角矩阵U。

（3）初始化迭代变量x_0。

（4）设置迭代次数N和容许误差ε。

（5）进行雅各比迭代计算，并输出每一步的迭代结果。

（6）判断迭代是否收敛，若收敛则输出最终结果，否则输出未收敛信息。

3. 运行实验程序，观察迭代过程和结果。

四、实验步骤1. 创建实验文件夹，打开MATLAB软件。

2. 编写实验程序，保存为“雅各比迭代法实验.m”。

3. 运行实验程序，观察迭代过程和结果。

4. 分析实验结果，验证雅各比迭代法的有效性和收敛性。

五、实验结果与分析1. 运行实验程序，得到以下迭代过程和结果：迭代次数 | 迭代结果---------|---------1 | x_1 = [0.3333, 0.3333]2 | x_2 = [0.3333, 0.3333]3 | x_3 = [0.3333, 0.3333]...N | x_N = [0.3333, 0.3333]2. 分析实验结果：（1）从实验结果可以看出，雅各比迭代法在求解线性方程组时，经过有限次迭代即可收敛。

数值分析迭代法数值分析实验报告 jacobi迭代和seidel迭代分析导读：就爱阅读网友为您分享以下“数值分析实验报告jacobi迭代和seidel迭代分析”资讯，希望对您有所帮助，感谢您对的支持!数值分析实验报告一、实验目的1、了解熟悉jacobi迭代法和seidel迭代法的解法2、将原理与matlab语言结合起来，编程解决问题3、分析实验结果二、实验题目设线性方程组为2 2 31 201 31 x1 1 7 x801240 2 x 103 10 15 x 1 4考察用jacobi迭代法和seidel迭代法求解该线性方程组的收敛情况。

如果收敛，给出误差满足x(k) x(k 1)10 4的解。

三、实验原理将A 作如下分解A D L U这里a11D a22 0 a ,L 21 anna n100 a n2 0 0 a12 a1n 0 0 0 a2n ,0 00 0 JacobiU迭代矩阵为BJ D 1(L U)Seidel迭代矩阵为BS (D L) 1U它们的迭代格式都可化为x(k 1) Bx(k) g则迭代格式对任何初值都瘦脸的充要条件是迭代矩阵谱半径 (B) max k 1,其中， k是矩阵B的n个特征值，k 1,2 ，n 1 k nJacobi迭代矩阵的特征方程为det( I B J) 0Seidel迭代矩阵的特征方程为det( I B S) det( (D L) U) 0四、实验内容用matlab编写计算jacobi迭代矩程序，建立m文件如下：function[M]=BJ(A)D=diag(diag(A));L=tril(-A)+D;U=triu(-A)+D;M=inv(D)*(L+U);输入：>> A=[20 1 -3 -1;3 18 0 7;-1 2 40 -2;1 0 -1 5];>> [M]=BJ(A)M:M =0 -0.0500 0.1500 0.0500-0.1667 0 0 -0.38890.0250 -0.0500 0 0.0500-0.2000 0 0.2000 00 0.050.150.05 0.166700 0.3889 则jacobi迭代矩阵为：BJ 0.025 0.0500.05 0.2 00.20用matlab求jacobi迭代矩阵的特征根的算法如下：>> A=[0 -0.05 0.15 -0.05;-0.67 0 0 -0.39;0.025 -0.05 0 0.05;-0.2 0 0.2 0];[V,D]=eig(A)V =-0.1892 + 0.0450i -0.1892 - 0.0450i -0.3812 -0.5005-0.9467 -0.9467 0.8867 0.5461 -0.1528 - 0.1181i -0.1528 + 0.1181i -0.2099 -0.0466 -0.1056 + 0.1325i -0.1056 - 0.1325i 0.1561 0.6701D =-0.1774 + 0.0864i 0 0 0 0 -0.1774 - 0.0864i 0 0 0 0 0.2194 0 0 0 0 0.1355 则最大特征根为：0.2194 则(BJ) 0.2194 1，所以jacobi迭代法收敛用matlab编程jacobi迭代法求根的算法：function [n,x]=jacobi(A,b,X,nm,w)%用雅克比迭代法求解方程组Ax=b%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w 为误差精度%输出：x为求得的方程组的解构成的列向量，n为迭代次数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)*(L+U); %计算迭代矩阵g=inv(D)*b; %计算迭代格式中的常数项%下面是迭代过程while n<=nmx=M*X+g; %用迭代格式进行迭代if norm(x-X,2)<wdisp(…迭代次数为‟);ndisp(…方程组的解为‟);xreturn;%上面：达到精度要求就结束程序，输出迭代次数和方程组的解endX=x;n=n+1;end%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）disp(…在最大迭代次数内不收敛!‟);输入数据：>> A=[20 1 -3 -1;3 18 0 7;-1 2 40 -2;1 0 -1 5];b=[1;2;10;-1];X=[0;0;0;0];>> nm=100;>> w=1e-4;>> [n,x]=jacobi1(A,b,X,nm,w)迭代次数为n =6方程组的解为x =0.06870.16450.2352-0.1667n =6x =0.06870.16450.2352-0.1667所以，满足精度的根为x = 0.0687 0.1645 0.2352用matlab编程计算seidel迭代矩阵算法为：function[M]=BS(A)D=diag(diag(A));L=tril(-A)+D;U=triu(-A)+D;M=inv(D-L)*U;输入：>> A=[20 1 -3 -1;3 18 0 7;-1 2 40 -2;1 0 -1 5]; >> [M]=BS(A)M =0 -0.0500 0.1500 0.0500 0 0.0083 -0.0250 -0.3972 -0.1667 迭代次数为6次0 -0.0017 0.0050 0.07110 0.0097 -0.0290 0.0042则seidel迭代矩阵为：0.150.05 0 0.05 00.0083 0.025 0.3972BS 0 0.00170.0050.0711 00.0097 0.0290.0042用matlab编程求得seidel矩阵的算法为：>> A=[0 -0.05 0.15 0.05;0 0.0083 -0.025 -0.3972;0 -0.0017 0.005 0.0711;0 0.0097 -0.029 0.0042];[V,D]=eig(A)V =1.0000 -0.1596 + 0.6750i -0.1596 - 0.6750i -0.91040 0.6965 0.6965 0.3924 0 -0.1250 + 0.0028i -0.1250 - 0.0028i 0.1313 0 0.0070 - 0.1348i 0.0070 + 0.1348i 0.0000D =0 0 0 0 0 0.0088 + 0.0768i 0 0 0 0 0.0088 - 0.0768i 0 0 0 0 -0.0001则特征根为 0 0.008+0.0768i 0.0088-0.0768i-0.0001则 (BS) 0.08 1，所以seidel迭代法收敛用seidel迭代法求根的算法为：function [n,x]=gaussseidel(A,b,X,nm,w)%用高斯-赛德尔迭代法求解方程组Ax=b%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w 为误差精度%输出：x为求得的方程组的解构成的列向量，n为迭代次数n=1;m=length(A);I=eye(m); %生成m*m阶的单位矩阵D=diag(diag(A)); %令A=D-L-U,计算矩阵DL=tril(-A)+D; %令A=D-L-U,计算矩阵LU=triu(-A)+D; %令A=D-L-U,计算矩阵U M=inv(D-L)*U; %计算迭代矩阵g=inv(I-inv(D)*L)*(inv(D)*b); %计算迭代格式中的常数项%下面是迭代过程while n<=nmx=M*X+g; %用迭代格式进行迭代if norm(x-X,2)<wdisp(…迭代次数为‟);ndisp(…方程组的解为‟);xreturn;%上面：达到精度要求就结束程序，输出迭代次数和方程组的解endX=x;n=n+1;end%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）disp(…在最大迭代次数内不收敛!‟);disp(…最大迭代次数后的结果为‟);x输入数据：>> A=[20 1 -3 -1;3 18 0 7;-1 2 40 -2;1 0 -1 5];b=[1;2;10;-1];X=[0;0;0;0];nm=100;w=1e-4;[n,x]=gaussseidel(A,b,X,nm,w) 迭代次数为n =5方程组的解为x =0.06870.16450.2352-0.1667n =5x =0.06870.16450.2352-0.166满足精度的根为x =0.0687 0.1645 0.2352 -0.1667 迭代次数为5次五、实验分析从实验过程可以看到求出满足精度的根，seidel迭代要比jacobi迭代快，这是因为在jacobi迭代计算中，计算向量x(k 1)时是按分量的角标由小到大依次计算的。

实验三一．实验题目分析用雅可比迭代法解线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡62525041-01-001-41-01-001-41-01-1-01-41-001-01-41-001-01-4654321x x x x x x 的收敛性，并求出使0001.02)()1(≤-+k k X X 得近似解及相应的迭代次数。

二．实验分析1.题目分析因为矩阵[]66⨯=ij a A 满足ij ni j j ii a a ≠=∑>1，所以该方程组有唯一解，雅可比迭代法收敛。

求使0001.02)()1(≤-+k k X X 的近似解，即求0001.0)(12)1(≤-∑=+ni k i k i x x 时的近似解。

2.代码分析 因为)6,2,1)((11⋅⋅⋅=-=∑≠=i x a b a x j n ij j ij i ii i 所以先给数组x[N]赋初值0，然后按公式计算，每一次算完x[N]后计算是否满足条件0001.0)(12)1(≤-∑=+n i k i k i x x， 若满足，则跳出循环，并输出近似解；反之，继续计算，直到迭代次数等于设定的允许的最大迭代次数。

三．代码实现#include<iostream>#include<cmath>using namespace std;#define N 6//方程阶数void main(){double a[N][N]={{4,-1,0,-1,0,0},{-1,4,-1,0,-1,0},{0,-1,4,-1,0,-1},{-1,0,-1,4,-1,0},{0,-1,0,-1,4,-1},{0,0,-1,0,-1,4}};double b[N]={0,5,-2,5,-2,6};double x[N]={0,0,0,0,0,0};double y[N];int k,i,j;int t=50;for(k=0;k<t;k++){ double e=0.0;for(i=0;i<N;i++){double m=0.0;for(j=0;j<N;j++){if(j!=i)m+=a[i][j]*x[j];}y[i]=(b[i]-m)/a[i][i];e+=pow(y[i]-x[i],2);cout<<"y["<<i+1<<"]="<<y[i]<<endl;}cout<<"e="<<e<<endl;cout<<"--------------------------"<<endl;if(sqrt(e)<=0.0001){cout<<"迭代次数k="<<k+1<<endl;for(i=0;i<N;i++){cout<<"y["<<i+1<<"]="<<y[i]<<endl;}break;}else{if(k==t-1)cout<<"计算失败！"<<endl;for(i=0;i<N;i++)x[i]=y[i];}}}四．运行结果五．实验小结 代码中需注意j n i j j ij x a m ∑≠==1与∑=+-=n i k i k i x x e 1)1()(的位置，m=每一行计算结束清零，e 每计算迭代完一次清零。

数学与计算科学学院《数值分析》课程设计题目：迭代法解线性方程组专业：信息与计算科学学号：*******-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 xk k =∞→，则*x 就是方程①或方程②的解。

此时迭代公式②是收敛的，否则称为发散的。

构造的迭代公式③是否收敛，取决于迭代矩阵B 的性1.雅可比迭代法基本原理设有方程组),,3,2,1(1n i b x aj 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 xk +==)()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 。

仿真平台与工具应用实践Jacobi迭代法求解线性方程组实验报告院系:专业班级:姓名:学号:指导老师:一、实验目的熟悉Jacobi迭代法原理；学习使用Jacobi迭代法求解线性方程组；编程实现该方法；二、实验内容应用Jacobi迭代法解如下线性方程组:, 要求计算精度为三、实验过程（1）、算法理论迭代格式的引出是依据迭代法的基本思想: 构造一个向量系列, 使其收敛至某个极限, 则就是要求的方程组的准确解。

Jacobi迭代将方程组:在假设, 改写成如果引用系数矩阵, 及向量, , ,方程组（1）和（2）分别可写为: 及, 这样就得到了迭代格式用迭代解方程组时, 就可任意取初值带入迭代可知式, 然后求。

但是, 比较大的时候, 写方程组和是很麻烦的, 如果直接由, 能直接得到, 就是矩阵与向量的运算了, 那么如何得到, 呢？实际上, 如果引进非奇异对角矩阵将分解成：要求的解, 实质上就有而是非奇异的, 所以存在, 从而有我们在这里不妨令就得到迭代格式：（2）算法框图（3）、算法程序m 文件:function x=jacobi(A,b,P,delta,n)N=length(b); %返回矩阵b的最大长度for k=1:nfor j=1:Nx(j)=(b(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(x'-P)); %求（x'-P）模的绝对值P=x';if(err<delta) %判断是否符合精度要求break;endendE=eye(N,N); %产生N行N列矩阵D=diag(diag(A));f=A*inv(D); %f是A乘D的逆矩阵B=E-f;Px=x';k,errBMATLAB代码:>> clear allA=[4, -1, 1;4, -8, 1;-2, 1, 5];b=[7, -21, 15]';P=[0,0,0]';x=jacobi(A,b,P,1e-7,20)（4）、算法实现用迭代法求解方程组:正常计算结果是2, 3, 4 , 下面是程序输出结果：P =2.00004.00003.0000k =17err =9.3859e-008B =0 -0.1250 -0.2000-1.0000 0 -0.20000.5000 0.1250 0x =2.00004.00003.0000四、实验体会五、MATLAB是非常实用的软件, 能够避免大量计算, 简化我们的工作, 带来便捷。

1.雅可比迭代源代码A=[6,-2,-1,-1;-2,18,-1,-1;-1,-1,6,-2;-1,-1,-1,12];b=[-16;6;9;43];x0=[0;0;0;0];it_max=500;eps=1e-6;[x,k]=jacobi(A,b,x0,eps,it_max)function [x,n]=jacobi(A,b,x0,eps,it_max)% 求线性方程组的Jacobi迭代法，调用格式为% [x, k] = jacobi(A,b,x0,eps,it_max)% 其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-6, % it_max 为最大迭代次数，默认为200% x 为线性方程组的解，k迭代次数if nargin ==3eps = 1.0e-6;M = 200;elseif nargin<3disp('输入参数数目不足3个');returnelseif nargin ==5M = it_max;endD = diag(diag(A));%求A的对角矩阵L = -tril(A,-1);%求A的下三角矩阵U = -triu(A,1);%求A的上三角矩阵B = D\(L+U);f = D\b;x = B*x0+f;n = 1;%迭代次数while norm(x-x0)>=epsx0 = x;x = B*x0+fn = n+1;if(n>=M)disp('Warning:迭代次数太多,可能不收敛!')return;endEnd2.高斯赛德尔迭代源代码A=[6,-2,-1,-1;-2,18,-1,-1;-1,-1,6,-2;-1,-1,-1,12];b=[-16;6;9;43];x0=[0;0;0;0];it_max=500;eps=1e-6;[x,k]=jacobi(A,b,x0,eps,it_max)function [x,n] = guaseidel(A,b,x0,eps,it_max)% 求线性方程组的Gauss-Seidel迭代法，调用格式为% [x, k] = guaseidel(A,b,x0,eps,it_max)% 其中, A 为线性方程组的系数矩阵，b 为常数项，eps 为精度要求，默认为1e-5, % it_max 为最大迭代次数，默认为100% x 为线性方程组的解，k迭代次数if nargin == 3eps = 1.0e-6it_max= 200elseif nargin == 4it_amx = 200elseif nargin <3disp('输入参数个数不足3个');return;endD = diag(diag(A));%求A的对角矩阵L = -tril(A,-1);%求A的下三角矩阵,不带对角线U = -triu(A,1);%求A的上三角矩阵G = (D-L)\U;f = (D-L)\b;x = G*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0 = x;x = G*x0+f;n = n+1;if(n>=it_max)disp('Warning:迭代次数太多,可能不收敛');return;endEnd3.方程组A=[6,-2,-1,-1;-2,18,-1,-1;-1,-1,6,-2;-1,-1,-1,12]; b=[-16;6;9;43];即4.MATLAB运行结果截图jacobi运行结果guaseidel运行结果jacobi迭代和guaseidel迭代从以上的运行结果比较来看，guaseidel的收敛性和收敛速度更快。

一、 实验目的与要求计算曲边梯形面积dx e I x⎰=20，可用内接梯形的面积n S 来进行逼近，用动画形象地描述这一逼近的过程（如果],[)(b a C x f ∈，按照定积分的定义，有∑⎰-=+=≈=10)()(n i n b a h ih a f S dx x f I ，其中n a b h -=。

这实际上是一个最简单的复化求积公式）实验要求：（1）用数值实验说明I S n n =∞→lim （其中2,,)(===b o a e x f x ）；（2）设计一副动画，模拟I S n →的过程一、 实验方案(程序源文件)a=input('请输入积分上限a= ');b=input('请输入积分上限b= ');n=[10,50,100,200,10000];m=length(n);for k=1:ms=0;i=1;h=(b-a)/n(k);g=exp(a);while i<n(k)y=exp(a+(i-1)*h);g=[g;y];s=s+y*h;i=i+1;endx=0:0.1:2;plot(x,exp(x),'r')hold onf=sym('2.71828^x');fprintf('曲边梯形面积')digits(5)I=int(f,'x',a,b)title(['y=exp(x),n=',num2str(n(k))])fprintf('内接梯形的面积')sx2=0.5*h:h:2-0.5*h;bar(x2,g,1.0,'c')pause(2)hold offend二、实验结果和数据处理(运行结果) 结果：请输入积分上限a= 0请输入积分上限b= 2曲边梯形面积I =6.3891内接梯形的面积s = 4.5615曲边梯形面积I = 6.3891内接梯形的面积s = 5.9782曲边梯形面积I =6.3891内接梯形的面积s = 6.1805曲边梯形面积I =6.3891内接梯形的面积s = 6.2840曲边梯形面积I =6.3891内接梯形的面积s = 6.3869曲边梯形面积I =6.3891内接梯形的面积s = 4.5615曲边梯形面积I =6.3891内接梯形的面积s = 5.9782曲边梯形面积I =6.3891内接梯形的面积s = 6.1805曲边梯形面积I =6.3891内接梯形的面积s = 6.2840曲边梯形面积I =6.3891内接梯形的面积s = 6.3869二、结论依图可知,当n趋于无穷大时，内接阶梯形面积趋近于曲边梯形的面积。

实验五矩阵的lu分解法，雅可比迭代法班级：学号：姓名：实验五矩阵的LU 分解法，雅可比迭代一、目的与要求：熟悉求解线性方程组的有关理论和方法；会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序； 通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。

二、实验内容：会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。

三、程序与实例 列主元高斯消去法算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+⨯表示 1) 消元过程 对k=1,2,…,n-1①选主元，找{}n ,,1k ,k i k Λ+∈使得k ,i k a =ik a ni k max≤≤②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。

③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，j i k j k a a ↔ j=k,┅,n+1④消元，对i=k+1, ┅,n 计算kk ik ik a a l /=对j=l+1, ┅,n+1计算kj ik ij ij a l a a -=2) 回代过程①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。

②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算ii ni j j ij n i i a x a a x /11,⎪⎪⎭⎫ ⎝⎛-=∑+=+程序与实例程序设计如下：#include <iostream> #include <cmath> using namespace std;void disp(double** p,int row,int col){ for(int i=0;i<row;i++){ for(int j=0;j<col;j++) cout<<p[i][j]<<' '; cout<<endl; } }void disp(double* q,int n){cout<<"====================================="<<endl; for(int i=0;i<n;i++)cout<<"X["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl; }void input(double** p,int row,int col){ for(int i=0;i<row;i++){ cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}int findMax(double** p,int start,int end){int max=start;for(int i=start;i<end;i++){if(abs(p[i][start])>abs(p[max][start]))max=i;}return max;}void swapRow(double** p,int one,int other,int col){double temp=0;for(int i=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}bool dispel(double** p,int row,int col){for(int i=0;i<row;i++){int flag=findMax(p,i,row); 用列主元消去法解方程⎪⎩⎪⎨⎧=++-=++-=++035.3643x .5072x .1835x .2137.2623x .43712x 347x .1 1.1833.555x 2.304x 0.101x 321321321例2 解方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧=++++=++++=++++=++++=++++-12.041.0F 1.02E 3.47D 1.04C 3.54B -6.301.0F2.01E 2.51D 4.04C 5.05B -8.531.0F 1.21E 2.92D 1.46C3.53B -20.071.0F 1.10E 4.48D 1.21C 4.93B -32.041.0F 1.55E 5.66D 2.40C 8.77B 计算结果如下B= C= D= E= F=矩阵直接三角分解法算法：将方程组A x=b 中的A 分解为A =LU ，其中L 为单位下三角矩阵，U 为上三角矩阵，则方程组A x=b 化为解2个方程组Ly =b ，Ux =y ，具体算法如下： ①对j=1,2,3,…,n 计算j j a u 11=对i=2,3,…,n 计算1111/a a l i i =②对k=1,2,3,…,n: a. 对j=k,k+1,…,n 计算∑-=-=11k q qj kq kj kj u l a ub. 对i=k+1,k+2,…,n 计算kk k q qk iq ik ik u u l a l /)(11∑-=-=③11b y =，对k=2,3,…,n 计算∑-=-=11k q q kq k k y l b y④nn n n u y x /=,对k=n-1,n-2,…,2,1计算kk nk q q kqk k u x uy x /)(1∑+=-=注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵[A ∣b ]=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++1,211,2222211,111211n n nn n n n n n n a a a a a a a a a a a a ΛM M M M M ΛΛ 施行算法②，③，此时U 的第n+1列元素即为y 。

实验五矩阵的lu分解法,雅可比迭代法班级:学号:姓名:实验五ﻩﻩ矩阵的L Ｕ分解法,雅可比迭代一、目的与要求:➢ 熟悉求解线性方程组的有关理论和方法;➢ 会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序; ➢ 通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。

二、实验内容:➢ 会编制列主元消去法、L Ｕ 分解法、雅可比及高斯—塞德尔迭代法德程序,进一步了解各种方法的优缺点。

三、程序与实例➢ 列主元高斯消去法算法：将方程用增广矩阵［Ａ∣b ］=（ij a )1n (n )+⨯表示1) 消元过程对k=1,2,…,n —１①选主元，找{}n ,,1k ,k i k +∈使得k ,i k a =ik a ni k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③.③如果k i k ≠,则交换第k 行与第k i 行对应元素位置,j i k j k a a ↔ j=k,┅，n+1④消元,对i=k+1, ┅,n 计算kk ik ik a a l /=对j=ｌ+1, ┅,n+1计算kj ik ij ij a l a a -=2) 回代过程①若0=nn a ，则矩阵Ａ奇异,程序结束;否则执行②。

②nn n n n a a x /1,+=;对i=n —1, ┅,2,1，计算ii n i j j ij n i i a x a a x /11,⎪⎪⎭⎫ ⎝⎛-=∑+=+程序与实例程序设计如下:#include ＜iｏstreaｍ>＃inｃlude <ｃmath>using ｎamespａce stｄ;voiｄ disｐ（ｄoｕｂle＊* p,ｉnt rｏｗ,inｔcｏl){for（ｉnt i=0;i<roｗ;i++)｛for(int j=0;ｊ〈ｃol；j＋＋)ｃout〈<p［i][j]〈〈' ';ｃｏut〈〈ｅndl;}｝ｖoid dｉsp（ｄｏuble＊q，ｉnt n){couｔ＜〈”==＝＝=＝==＝===＝=＝=＝===＝===============＝”〈<endl;for(int i=0;ｉ〈n;ｉ++）coｕt〈<”X［”〈<ｉ+１<<"]=”＜<q[ｉ］〈〈endl；ｃouｔ〈〈"===＝=＝＝===＝======＝＝＝=＝==============＝”<<eｎdl；｝ｖoid iｎｐut(doｕble＊* p,int ｒow,int coｌ){fｏr（int i=０;i〈row;i＋+){cｏut＜<”输入第"<<ｉ＋１〈〈”行：";foｒ(int j＝0;j<ｃol；ｊ++)ｃin＞>p[i][j］;}}ｉnt findＭax(double＊＊ p，int stａｒt,inｔ enｄ){ｉnt max=starｔ；for（iｎt i=sｔart;ｉ＜ｅnd;i++){if（abs(p[i]［ｓｔarｔ])>aｂs(p［max］［stａrt］))ｍax=i;｝reｔuｒn max；｝void ｓwapRow（douｂlｅ** p，int one,int other,int col){ｄouｂle ｔｅmp=0;ｆor(int i＝0；i〈ｃol;i+＋){temp＝p[onｅ］[ｉ］;p［ｏne]［i］＝p[ｏtheｒ］［i];p[otheｒ][ｉ]=ｔemp；}}booｌｄispel（doｕble*＊ p,ｉnt row，iｎt cｏｌ){ｆｏr(ｉnt i=０；i〈ｒow;ｉ++)｛ｉnt flａｇ=findＭａｘ(ｐ,i,ｒoｗ）；/／找列主元行号ｉf(p[ｆlag]［ｉ]==0）rｅｔurn fａｌse；swａpRoｗ（p，i,ｆlaｇ,ｃｏl)；／/交换行for（int j=i＋1;j〈row；ｊ++)｛ｄouｂle eｌem＝ｐ[j][ｉ］/p［i][i]; ／／消元因子p［j]［i]=0；ｆor(ｉnt ｋ=i＋1;k<cｏl；ｋ++）{ｐ[ｊ］[k]—＝（eleｍ*p[i］[ｋ］）；｝｝}ｒｅtuｒn tｒｕe；｝douｂle sｕmRｏw(ｄouｂlｅ＊＊ p,douｂle＊ q，int ｒoｗ，int coｌ){ doubｌe suｍ=０；ｆoｒ(ｉnt i=0;i＜col-1;i++）{ｉf(ｉ==row)ｃontinuｅ;ｓum+=(q［i］＊p[rｏw][i］);}retuｒn ｓum;}ｖｏid back(double*＊ p,ｉnｔ rｏw，ｉｎt col，double＊ｑ）｛ for（ｉnt i=roｗ—１;ｉ＞＝0;i－—){ｑ[i]=(ｐ[ｉ]［col—1］－suｍRow(p，q，i，ｃol））/p[i］［i]；｝}inｔ main（)｛cout〈<"Inpｕt ｎ:";inｔ n; //方阵的大小cin〉>ｎ;ｄoubｌe **p=ｎew ｄoｕbｌｅ* ［n]；for（inｔ i＝0;i〈n;i＋+）{p[i]=neｗ double [ｎ+1］;｝in ｐut(p ，n,n+1);i ｆ(!ｄispel(p,n,n ＋１））{co ｕt<〈”奇异＂<<en ｄｌ;ｒeturn 0;}ｄoub ｌｅ* ｑ=ｎew doubl ｅ［n];fo ｒ(int ｉ=0;i ＜n ；ｉ++）ｑ[ｉ]＝0;back （p,n,n+1，ｑ）;disp(q,n)；delete[] q;ｆｏr （i ｎt i=0;i ＜ｎ；ｉ＋＋)del ｅte[] ｐ[ｉ];d ｅlete ［］ p;｝１。

实验五矩阵的lu分解法，雅可比迭代法班级：学号：：实验五矩阵的LU分解法，雅可比迭代一、目的与要求：➢熟悉求解线性方程组的有关理论和方法；➢会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序；➢通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。

二、实验容：➢会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。

三、程序与实例➢列主元高斯消去法算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+⨯表示1) 消元过程对k=1,2,…,n-1①选主元，找{}n ,,1k ,k i k +∈使得k ,i k a =ik a ni k max ≤≤ ②如果0a k ,i k =，则矩阵A 奇异，程序结束；否则执行③。

③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，j i k j k a a ↔ j=k,┅,n+1④消元，对i=k+1, ┅,n 计算kk ik ik a a l /=对j=l+1, ┅,n+1计算kj ik ij ij a l a a -=2) 回代过程①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。

②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算ii ni j j ij n i i a x a a x /11,⎪⎪⎭⎫ ⎝⎛-=∑+=+程序与实例程序设计如下： #include <iostream>#include <cmath>using namespace std;void disp(double** p,int row,int col){for(int i=0;i<row;i++){for(int j=0;j<col;j++)cout<<p[i][j]<<' ';cout<<endl;}}void disp(double* q,int n){cout<<"====================================="<<endl;for(int i=0;i<n;i++)cout<<"X["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl; }void input(double** p,int row,int col){for(int i=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(int j=0;j<col;j++)cin>>p[i][j];}}int findMax(double** p,int start,int end){int max=start;for(int i=start;i<end;i++){if(abs(p[i][start])>abs(p[max][start]))max=i;}return max;}void swapRow(double** p,int one,int other,int col){ double temp=0;for(int i=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}bool dispel(double** p,int row,int col){for(int i=0;i<row;i++){int flag=findMax(p,i,row); //找列主元行号if(p[flag][i]==0) return false;swapRow(p,i,flag,col); //交换行for(int j=i+1;j<row;j++){double elem=p[j][i]/p[i][i]; //消元因子p[j][i]=0;for(int k=i+1;k<col;k++){p[j][k]-=(elem*p[i][k]);}}}return true;}double sumRow(double** p,double* q,int row,int col){ double sum=0;for(int i=0;i<col-1;i++){if(i==row)continue;sum+=(q[i]*p[row][i]);}return sum;}void back(double** p,int row,int col,double* q){ for(int i=row-1;i>=0;i--){q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i];}}int main(){cout<<"Input n:";int n; //方阵的大小cin>>n;double **p=new double* [n];for(int i=0;i<n;i++){p[i]=new double [n+1];}input(p,n,n+1);if(!dispel(p,n,n+1)){cout<<"奇异"<<endl;return 0;}double* q=new double[n];for(int i=0;i<n;i++)q[i]=0;back(p,n,n+1,q);disp(q,n);delete[] q;for(int i=0;i<n;i++)delete[] p[i];delete[] p;}1. 用列主元消去法解方程⎪⎩⎪⎨⎧=++-=++-=++035.3643x .5072x .1835x .2137.2623x .43712x 347x .1 1.1833.555x 2.304x 0.101x 321321321例2 解方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧=++++=++++=++++=++++=++++-12.041.0F 1.02E 3.47D 1.04C 3.54B -6.301.0F2.01E 2.51D 4.04C 5.05B -8.531.0F 1.21E2.92D 1.46C3.53B -20.071.0F 1.10E4.48D 1.21C 4.93B -32.041.0F 1.55E5.66D 2.40C 8.77B 计算结果如下B=-1.461954C= 1.458125D=-6.004824 E=-2.209018F= 14.719421➢ 矩阵直接三角分解法算法：将方程组A x=b 中的A 分解为A =LU ，其中L 为单位下三角矩阵，U 为上三角矩阵，则方程组A x=b 化为解2个方程组Ly =b ，Ux =y ，具体算法如下：①对j=1,2,3,…,n 计算j j a u 11=对i=2,3,…,n 计算1111/a a l i i =②对k=1,2,3,…,n:a. 对j=k,k+1,…,n 计算∑-=-=11k q qj kq kj kj u l a ub. 对i=k+1,k+2,…,n 计算kk k q qk iq ik ik u u l a l /)(11∑-=-=③11b y =，对k=2,3,…,n 计算∑-=-=11k q q kq k k y l b y④nn n n u y x /=,对k=n-1,n-2,…,2,1计算kk n k q q kq k k u x u y x /)(1∑+=-=注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵[A ∣b ]=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+++1,211,2222211,111211n n nn n n n nn n a a a a a a a a a a a a 施行算法②，③，此时U 的第n+1列元素即为y 。