大连理工大学矩阵大作业
大连理工大学《矩阵与数值分析》2005-2009年真题答案
大 连 理 工 大 学课 程 名 称： 计算方法 试卷： A 考试形式： 闭卷 授课院（系）： 数学系 考试日期： 2005 年 12 月 12 日 试卷共 7 页一二三四五 六 七 总分 标准分 得 分装 一、填空（共30分，每空1.5分）（1）误差的来源主要有 、 、 、 .（2）要使 7459666.760=的近似值a 的相对误差限不超过310-，应至少取 位有效数字, 此时的近似值a = .订 （3）设⎪⎪⎭⎫⎝⎛--=4224A , 则1A ＝ , 2A ＝ , ∞A ＝ , F A ＝ ,谱半径)(A ρ＝ , 2-条件数)(2A cond ＝ , 奇异值为 .线 （4）设44⨯∈CA ,特征值3,24321====λλλλ,特征值2是半单的，而特征值3是亏损的，则A 的Jordan 标准型=J.（5）已知x x x f 3)(2-=，则=-]1,0,1[f ，=-]3,1,0,1[f .（6）求01)(3=-+=x x x f 在5.0=x 附近的根α的Newton 迭代公式是：，其收敛阶 . （7）计算u u 5-=')10(≤≤t , 1)0(=u 的数值解的Euler 求解公式为 . 为使计算保持绝对稳定性, 步长h 的取值范围 .二、（12分）求矩阵⎪⎪⎪⎭⎫ ⎝⎛=820251014A 的Doolittle 分解和Cholesky 分解，并求解⎪⎪⎪⎭⎫ ⎝⎛=1085Ax .三、（6分）求矩阵⎪⎪⎪⎭⎫ ⎝⎛=622292221A 的QR 分解（Q 可表示为两个矩阵的乘积）.四、（12分）根据迭代法f Bx x k k +=+)()1(对任意)0(x 和f 均收敛的充要条件为1)(<B ρ, 证明若线性方程组b Ax =中的A 为严格对角占优矩阵, 则Jacobi 法和G-S 法均收敛.五、（12分）求满足下列插值条件的分段三次多项式（]0,3[-和]1,0[）, 并验证它是不是三次样条函数.27)3(-=-f , 8)2(-=-f , 1)1(-=-f , 0)0(=f , ]0,3[-∈x ；0)0(=f , 0)0(='f , 0)1(=f , 1)1(='f , ]1,0[∈x .六、（10分）证明线性二步法])13()3[(4)1(212n n n n n f b f b hbu u b u +++=--++++, 当1-≠b 时为二阶方法,1-=b 时为三阶方法, 并给出1-=b 时的局部截断误差主项.七、（18分）求]1,1[-上以1)(≡x ρ为权函数的标准正交多项式系)(0x ψ, )(1x ψ, )(2x ψ, 并由此求3x ])1,1[(-∈x 的二次最佳平方逼近多项式, 构造Gauss 型求积公式⎰-+≈111100)()()(x f A x f A dx x f , 并验证其代数精度.大 连 理 工 大 学课 程 名 称： 计算方法 试卷： A 考试形式： 闭卷 授课院（系）： 数学系 考试日期： 2006 年 12 月 11 日 试卷共 8 页一二三四五 六 七 八 总分 标准分 得 分装订 一、填空（共30分，每空2分）线 （1）误差的来源主要有 .（2）按四舍五入的原则，取 69041575.422= 具有四位有效数字的近似值 a = ，则绝对误差界为 ，相对误差界为 .（3）矩阵算子范数M A ||||和谱半径)(A ρ的关系为： ，和 .（4）设44⨯∈CA ,特征值3,24321====λλλλ,特征值2是半单的，而特征值3是亏损的，则A 的Jordan 标准型=J.（5）已知x x x f 3)(2-=，则=]1,0[f ，=-]1,0,1[f .（6）求01)(3=-+=x x x f 在5.0=x 附近的根α的Newton 迭代公式是：.（7）使用Aitken 加速迭代格式)(1-=k k x x ϕ得到的Steffensen 迭代格式为：，对幂法数列}{k m 的加速公式为：.（8）1+n 点的Newton-Cotes 求积公式∑==nk k k n x f A f I 0)()(的最高代数精度为.（9）计算u u 7-=')10(≤≤t , 1)0(=u 的数值解的Euler 求解公式为 ，为使计算保持绝对稳定性, 步长h 的取值范围 .二、（10分） 设⎪⎪⎭⎫ ⎝⎛--=4224A , 计算1A ,2A ,∞A ,F A , 谱半径)(A ρ, 2-条件数)(2A cond , 和奇异值.三、（10分）求矩阵⎪⎪⎪⎭⎫ ⎝⎛=820251014A 的Doolittle 分解和Cholesky 分解.四、（4分）求Householder 变换矩阵将向量⎪⎪⎪⎭⎫ ⎝⎛=221x 化为向量⎪⎪⎪⎭⎫ ⎝⎛=003y .五、（12分）写出解线性方程组的Jacobi 法，G-S 法和超松弛（SOR ）法的矩阵表示形式，并根据迭代法f Bx x k k +=+)()1(对任意)0(x 和f 均收敛的充要条件为1)(<B ρ, 证明若线性方程组b Ax =中的A 为严格对角占优矩阵, 则超松弛（SOR ）法当松弛因子]1,0(∈ω时收敛.六、（12分）求满足下列插值条件的分段三次多项式（]0,3[-和]1,0[）, 并验证它是不是三次样条函数. 27)3(-=-f , 8)2(-=-f , 1)1(-=-f , 0)0(=f , ]0,3[-∈x ；0)0(=f , 0)0(='f , 0)1(=f , 1)1(='f , ]1,0[∈x .七、（12分）证明区间],[b a 上关于权函数)(x ρ的Gauss 型求积公式∑==nk k k n x f A f I 0)()(中的系数⎰=bak k dx x l x A )()(ρ，其中)(x l k 为关于求积节点n x x x ,,10的n 次Lagrange 插值基函数，n k ,1,0=. 另求]1,1[-上以1)(≡x ρ为权函数的二次正交多项式)(2x ψ, 并由此构造Gauss型求积公式⎰-+≈111100)()()(x f A x f A dx x f .八、（10分）证明线性二步法])13()3[(4)1(212n n n n n f b f b hbu u b u +++=--++++, 当1-≠b 时为二阶方法, 1-=b 时为三阶方法, 并给出1-=b 时的局部截断误差主项.大连理工大学应用数学系数学与应用数学专业2005级试A 卷答案课 程 名 称： 计算方法 授课院 (系)： 应 用 数 学 系 考 试 日 期：2007年11 月 日 试卷共 6 页一 二 三 四 五 六 七 八 九 十 总分标准分 42 8 15 15 15 5 / / / / 100 得 分一、填空（每一空2分，共42分）1．为了减少运算次数，应将表达式.543242161718141311681x x x x x x x x -+---++- 改写为()()()()()()()1816011314181716-+++---+-x x x x x x x x x ；2．给定3个求积节点：00=x ，5.01=x 和12=x ，则用复化梯形公式计算积分dxe x ⎰-12求得的近似值为()15.02141--++e e ， 用Simpson 公式求得的近似值为()15.04161--++e e 。
大连理工大学《矩阵与数值分析》学习指导与课后参考答案第三章、逐次逼近法
第三章 逐次逼近法1.1内容提要1、一元迭代法x n+1=φ(x n )收敛条件为：1）映内性x ∈[a,b]，φ(x) ∈[a,b] 2）压缩性∣φ(x) -φ(y)∣≤L ∣x-y ∣其中L ＜1，此时φ为压缩算子，在不断的迭代中，就可以得到最终的不动点集。
由微分中值定理，如果∣φ’∣≤L ＜1，显然它一定满足压缩性条件。
2、多元迭代法x n+1=φ(x n )收敛条件为：1）映内性x n ∈Ω，φ(x n ) ∈Ω 2）压缩性ρ（▽φ）＜1，其中▽φ为x n 处的梯度矩阵，此时φ为压缩算子，在不断的迭代中，就可以得到最终的不动点集。
3、当φ(x )= Bx+f 时，收敛条件为，ρ（B ）＜1，此时x n+1= Bx n +f ，在不断的迭代中，就可以得到线性方程组的解。
4、线性方程组的迭代解法，先作矩阵变换 U L D A --= Jacobi 迭代公式的矩阵形式 f Bx b D x U L D x n n n +=++=--+111)(Gauss-Seidel 迭代公式的矩阵形式 f Bx b L D Ux L D x n n n +=-+-=--+111)()( 超松弛迭代法公式的矩阵形式f Bx b L D x U D L D x k k k +=-++--=--+ωωωωω111)(])1[()(三种迭代方法当1)(<B ρ时都收敛。
5、线性方程组的迭代解法，如果A 严格对角占优，则Jacob 法和Gauss-Seidel 法都收敛。
6、线性方程组的迭代解法，如果A 不可约对角占优，则Gauss-Seidel 法收敛。
7、Newton 迭代法，单根为二阶收敛 2211'''21lim)(2)(lim---∞→+∞→--=-==--k k k k k k k k x x x x f f c x x ξξαα8、Newton 法迭代时，遇到重根，迭代变成线性收敛，如果知道重数m ， )()('1k k k k x f x f m x x -=+仍为二阶收敛 9、弦割法)()())((111--+---=k k k k k k k x f x f x x x f x x 的收敛阶为1.618，分半法的收敛速度为（b-a ）/2n-110、Aitken 加速公式11211112)(),(),(+----+-+--+---+---===k k k k k k k k k k k x x x x x x x x x x x ϕϕ1.2 典型例题分析1、证明如果A 严格对角占优，则Jacob 法和Gauss-Seidel 法都收敛。
大连理工大学矩阵与数值分析上机作业代码
T
方程组，并比较计算结果。
1.1 程序
（1）高斯消元法 n=10 时， >> [A,b]=CreateMatrix(10) A= 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6 8 0 0 0 0 0 0 0 0 1 6
1.3 M 文件
.m 1.3.1 CreateMatrix CreateMatrix.m function [A,b]=CreateMatrix(n) %用于存放习题1的题目信息，并建立构造题目中矩阵的函数 %对矩阵A赋值 A1=6*ones(1,n); A2=ones(1,n-1); A3=8*ones(1,n-1); A=diag(A1)+diag(A2,1)+diag(A3,-1); %对向量b赋值 b=15*ones(n,1); b(1)=7; b(n)=14;
10
，迭代次数上限取默认值
50，使用 Jacobi 法进行迭代。 >> test2 >> b=ones(20,1) >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2227 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2227 0.2223 0.2159 0.2327 0.2766
大连理工矩阵上机作业
第一题Lagrange插值函数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;endx0=[1:10];y0=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76 .777];lagrange(x0,y0,17)ans=-1.9516e+12x=[1:0.1:10];x=x';plot(x0,y0,'r');hold onplot(x,y,'k');legend('原函数','拟合函数')拟合图像如下拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项式会因为误差的不断累积，导致龙格现象的发生。
第二题function fun =nihe(n)m=[67.052*10^6,68.008*10^6,69.803*10^6,72.024*10^6,73.400*10^6,72.063 *10^6,74.669*10^6,74.487*10^6,74.065*10^6,76.777*10^6];w=[1,2,3,4,5,6,7,8,9,10];d1=0;d2=0;d3=0;y1=polyfit(m,w,1);y2=polyfit(m,w,2);y3=polyfit(m,w,3);y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);f1=subs(y1,17);f2=subs(y2,17);f3=subs(y3,17);for p=1:10;d1=d1+(subs(y1,w(p))-m(p))^2;d2=d2+(subs(y2,w(p))-m(p))^2;d3=d3+(subs(y3,w(p))-m(p))^2;endd1=sqrt(d1);d2=sqrt(d2);d3=sqrt(d3);fun=[f1 f2 f3;d2 d3 d4];return;结果三次函数的均方误差最小，拟合的最好。
数据结构作业答案(大连理工大学)
作业1. 线性表编程作业：1．将顺序表逆置，要求用最少的附加空间。
参考答案#include <>#include <>#include <>#define LIST_INIT_SIZE 100#define LISTINCREMENT 10#define TRUE 1#define FALSE 0#define OK 1#define ERROR 0#define INFEASIBLE -1#define OVERFLOW -2typedef int Status;typedef int ElemType;typedef struct{ ElemType *elem;int length;int listsize;}SqList;立单链表 ");printf("2.取元素值 ");printf("3.查找 \n");printf("4.插入 ");printf("5.删除 ");printf("6.显示\n");printf("7.删除大于mink且小于maxk的元素值 ");printf("8.就地升序排序\n");printf("9.就地逆置 ");printf("a.有序表插入 ");printf("q.退出\n");printf("\n请选择操作：");fflush(stdin);scanf("%c",&choice);switch(choice){case '1': printf("请输入单链表中结点个数：");scanf("%d",&n);Create_L2(L,n);break;case '2': printf("请输入元素位序：");scanf("%d",&i);GetElem_L(L,i,e);printf("元素值为：%d\n",e);break;case '3': printf("请输入要查找的元素：");scanf("%d",&e);if(dlbcz(L,e))printf("查找成功！");elseprintf("查找失败。
大连理工大学矩阵第二章(矩阵变换和计算)
第二章 矩阵变换和计算一、内容提要本章以矩阵的各种分解变换为主要内容，介绍数值线性代数中的两个基本问题：线性方程组的求解和特征系统的计算，属于算法中的直接法。
基本思想为将计算复杂的一般矩阵分解为较容易计算的三角形矩阵. 要求掌握Gauss （列主元）消去法、矩阵的（带列主元的）LU 分解、平方根法、追赶法、条件数与误差分析、QR 分解、Shur 分解、Jordan 分解和奇异值分解.(一) 矩阵的三角分解及其应用 1．矩阵的三角分解及其应用考虑一个n 阶线性方程组b Ax =的求解，当系数矩阵具有如下三种特殊形状：对角矩阵D ，下三角矩阵L 和上三角矩阵U ，这时方程的求解将会变得简单. ⎪⎪⎪⎪⎪⎭⎫⎝⎛=n d dd D21, ⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=nnn n l l l l l l L21222111, ⎪⎪⎪⎪⎪⎭⎫⎝⎛=nn n n u u u u u u U22212111. 对于b Dx =，可得解为i i i d b x /=,n i ,,2,1 =.对于b Lx =，可得解为1111/l b x =,ii i k k iki i l x lb x /)(11∑-=-=,n i ,,3,2 =.对于b Ux =，可得解为nn n n l b x /=,ii ni k k iki i l x lb x /)(1∑+=-=,1,,2,1 --=n n i .虽然对角矩阵的计算最为简单，但是过于特殊，任意非奇异矩阵并不都能对角化，因此较为普适的方法是对矩阵进行三角分解.1）．Gauss 消去法只通过一系列的初等行变换将增广矩阵)|(b A 化成上三角矩阵)|(c U ，然后通过回代求与b Ax =同解的上三角方程组c Ux =的解．其中第k 步消元过程中，在第1-k 步得到的矩阵)1(-k A 的主对角元素)1(-k kka 称为主元.从)1(-k A 的第j 行减去第k 行的倍数)1()1(--=k kkk jkjk a a l （n j k ≤<）称为行乘数（子）.2）．矩阵A 的LU 分解对于n 阶方阵A ，如果存在n 阶单位下三角矩阵L 和n 阶上三角矩阵U ，使得LU A =, 则称其为矩阵A 的LU 分解，也称为Doolittle 分解．Gauss 消去法对应的矩阵形式即为LU 分解, 其中L 为所有行乘子组成的单位下三角矩阵, U 为Gauss 消去法结束后得到的上三角矩阵. 原方程组b Ax =分解为两个三角形方程组⎩⎨⎧==yUx b Ly .3）．矩阵LU 分解的的存在和唯一性如果n 阶矩阵A 的各阶顺序主子式),,2,1(n k k =D 均不为零, 则必有单位下三角矩阵L 和上三角矩阵U ，使得LU A =, 而且L 和U 是唯一存在的．4）．Gauss 列主元消去法矩阵每一列主对角元以下（含主对角元）的元素中, 绝对值最大的数称为列主元. 为避免小主元作除数、或0作分母，在消元过程中，每一步都按列选主元的Guass 消去法称为Gauss 列主元消去法．由于选取列主元使得每一个行乘子均为模不超过1的数，因此它避免了出现大的行乘子而引起的有效数字的损失.5）．带列主元的LU 分解Gauss 列主元消去法对应的矩阵形式即为带列主元的LU 分解,选主元的过程即为矩阵的行置换. 因此, 对任意n 阶矩阵A ，均存在置换矩阵P 、单位下三角矩阵L 和上三角矩阵U ，使得LU PA =．由于选列主元的方式不唯一, 因此置换矩阵P 也是不唯一的. 原方程组b Ax =两边同时乘以矩阵P 得到Pb PAx =, 再分解为两个三角形方程组⎩⎨⎧==y Ux PbLy .5）．平方根法(对称矩阵的Cholesky 分解)对任意n 阶对称正定矩阵A ，均存在下三角矩阵L 使T LL A =，称其为对称正定矩阵A 的Cholesky 分解. 进一步地, 如果规定L 的对角元为正数，则L 是唯一确定的．原方程组b Ax =分解为两个三角形方程组⎩⎨⎧==y x L bLy T .利用矩阵乘法规则和L 的下三角结构可得21112⎪⎪⎭⎫ ⎝⎛-=∑-=j k jkjj jjla l , jj j k jk ik ij ij l l l a l /11⎪⎪⎭⎫⎝⎛-=∑-=, i=j +1, j +2,…,n , j =1,2,…,n . 计算次序为nn n n l l l l l l l ,,,,,,,,,2322212111 ．由于jj jk a l ≤，k =1,2,…,j ．因此在分解过程中L 的元素的数量级不会增长，故平方根法通常是数值稳定的，不必选主元．6）．求解三对角矩阵的追赶法 对于三对角矩阵⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=---n nn n n b a c b a c b a c b 11122211A , 它的LU 分解可以得到两个只有两条对角元素非零的三角形矩阵 ⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=--n n n nu d u d u d u l l l 11221132,1111U L . 其中⎪⎪⎩⎪⎪⎨⎧=-====-==--n i c l b u n i u a l b u n i c d i i i i i i i i i ,,3,2,,,3,2,/1,,2,1,1111计算次序是n n u l u l u l u →→→→→→→ 33221. 原方程组b Ax =分解为两个三角形方程组⎩⎨⎧==y Ux b Ly . 计算公式为n i y l b y b y i i i i ,,3,2,,111 =-==-,.1,,2,1,/)(,/1 --=-==+n n i u x c y x u y x i i i i i nn n该计算公式称为求解三对角形方程组的追赶法．当A 严格对角占优时，方程组b Ax =可用追赶法求解, 解存在唯一且数值稳定．7）．矩阵的条件数设A 为非奇异矩阵，⋅为矩阵的算子范数，称1)(cond -=A A A 为矩阵A 的条件数．矩阵的条件数是线性方程组b Ax =, 当A 或b 的元素发生微小变化，引起方程组解的变化的定量描述, 因此是刻画矩阵和方程组性态的量. 条件数越大, 矩阵和方程组越为病态, 反之越小为良态.常用的矩阵条件数为∞-条件数: ∞-∞∞=1)(cond AA A ，1-条件数: 1111)(cond -=AAA ，2-条件数: )()()(cond mi n max 2122A A A A AAA HHλλ==-．矩阵的条件数具有如下的性质： (1) 1)(cond ≥A ;(2) )(cond )(cond 1-=A A ;(3) )(cond )(cond A A =α，0≠α，R ∈α;(4) 如果U 为正交矩阵，则1)(cond 2=U ，)(cond )(cond )(cond 222A AU UA ==.一般情况下，系数矩阵和右端项的扰动对解的影响为定理 2.5 设b Ax =，A 为非奇异矩阵，b 为非零向量且A 和b 均有扰动．若A 的扰动δA 非常小,使得11<-A A δ，则)()(cond 1)(cond bδb AδA AA A A xδx +-≤δ.关于近似解的余量与它的相对误差间的关系有定理2.6 设b Ax =，A 为非奇异矩阵，b 为非零向量，则方程组近似解x ~的事后估计式为bx A b A xx x bx A b A ~)cond(~~)cond(1-≤-≤-.其中称x A b ~-为近似解x ~的余量，简称余量。
大连理工大学矩阵与数值分析上机作业
end
case2%2-范数
fori=1:n
s=s+x(i)^2;
end
s=sqrt(s);
caseinf%无穷-范数
s=max(abs(x));
end
计算向量x，y的范数
Test1.m
clearall;
clc;
n1=10;n2=100;n3=1000;
x1=1./[1:n1]';x2=1./[1:n2]';x3=1./[1:n3]';
xlabel('x');ylabel('p(x)');
运行结果：
x=2的邻域：
x =
1.6000 1.8000 2.0000 2.2000 2.4000
相应多项式p值：
p =
1.0e-003 *
-0.2621 -0.0005 0 0.0005 0.2621
p(x)在 [1.95,20.5]上的图像
程序：
[L,U]=LUDe.(A);%LU分解
xLU=U\(L\b)
disp('利用PLU分解方程组的解:');
[P,L,U] =PLUDe.(A);%PLU分解
xPLU=U\(L\(P\b))
%求解A的逆矩阵
disp('A的准确逆矩阵:');
InvA=inv(A)
InvAL=zeros(n);%利用LU分解求A的逆矩阵
0 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625
0 0 0 0.5000 -0.2500 -0.1250 -0.1250
0 0 0 0 0.5000 -0.2500 -0.2500
大连理工大学矩阵分析matlab上机作业
x(i)=1/i; %按要求给向量 x 赋值，其值递减 end normx1=norm(x,1); %求解向量 x 的 1 范数 normx1 normx2=norm(x,2); %求解向量 x 的 2 范数 normx2 normxinf=norm(x,inf); %求解向量 x 的无穷范数 normxinf normy1=norm(y,1); %求解向量 y 的 1 范数 normy1 normy2=norm(y,2); %求解向量 y 的 2 范数 normy2 normyinf=norm(y,inf); %求解向量 y 的无穷范数 normyinf z1=[normx1,normx2,normxinf]; z2=[normy1,normy2,normyinf]; end
for i=2:n
for j=i:n U(i,j)=A(i,j)-L(i,1:i-1)*U(1:i-1,j);
式
%Doolittle 分解计算上三角矩阵的公
L(j,i)=(A(j,i)-L(j,1:i-1)*U(1:i-1,i))/U(i,i); %Doolittle 分解计算下三角矩 阵的公式
end
1 1 1 ������ x = (1, 2 , 3 , … , ������) ,
������ = (1,2, … , ������)������.
对n = 10，100，1000甚至更大的n计算其范数，你会发现什么结果?你能否修改
你的程序使得计算结果相对精确呢？
1.1 源代码
function [z1,z2]=norm_vector(n) %向量 z1 的值为向量 x 的是三种范数，向量 z2 的值为向量 y 的三 种范数，n 为输入参数
大连理工大学矩阵与数值分析试卷-2013
1 0 0 0
3 ⎞ ⎛2 5 ⎟ T ⎟ ; LL 分解中 L= ⎜ ⎜3 4 ⎜ − ⎟ ⎟ ⎝2 5⎠
1 1 2 2
0 ⎞ ⎟ 7 ⎟。 ⎟ 2 ⎠
Gauss 求 积 公 式 ， 则
1 ∫ x + 1 f (x ) dx ≈ A f (x ) + A f (x ) + A f (x ) 为
2）为使二点数值求积公式 积节点和求积系数应为 (A) x0 = −
∫
1
f ( x) 1 − x2
.
−1
dx ≈ A0 f ( x0 ) + A1 f ( x1 ) 具有最高的代数精度，其求
B
2 2 π 1 1 1 , x1 = ; A0 = A1 = ； (B) x0 = − , x1 = ; A0 = A1 = ； 2 2 2 2 2 2
⎛ ⎜ 即 V = ( v1 v2 ) = ⎜ ⎜ ⎜ ⎝ ⎛ ⎜ V1 = V = ⎜ ⎜ ⎜ ⎝ 1 2 1 2 1 ⎞ ⎛ ⎟ ⎜ 2⎟ 或 V = ( v1 v2 ) = ⎜ −1 ⎟ ⎜ ⎟ ⎜ 2⎠ ⎝ 1 2 1 2 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ −1 ⎞ 2⎟ ⎟ ，因 rank(A)=1，故有 1 ⎟ ⎟ 2⎠ 1 ⎞ ⎛ 1 ⎞ ⎟ ⎜ ⎟ 2⎟ (1) = ⎜ 2 ⎟ , 由 U = (U1U 2 ) , 则 1 ⎟ ⎜ 1 ⎟ ⎟ ⎜ ⎟ 2⎠ ⎝ 2⎠
17）. 为了减少运算次数，应将表达式.
4 x3 − 3x 2 − 2 x − 1 改写为 x4 + x2 + x − 1
( ( 4 x − 3) x − 2 ) x − 1 ； ( ( ( x + 0 ) x + 1) x + 1) x − 1
大连理工_2012矩阵与数值分析大作业
矩阵与数值分析学生：学号：任课老师：金光日教学班号：（2）班院系：电子信息与电气工程学部《矩阵与数值分析》课程数值实验题目1.给定n 阶方程组A x b =，其中6186186186A ⎛⎫ ⎪ ⎪⎪= ⎪ ⎪ ⎪⎝⎭，7151514b ⎛⎫ ⎪⎪ ⎪= ⎪ ⎪⎪⎝⎭则方程组有解(1,1,,1)T x = 。
对10n =和84n =，分别用Gauss 消去法和列主元消去法解方程组，并比较计算结果。
1答： 程序1. Gauss 消元法function x=DelGauss(A,b) % Gauss 消去法 [n,m]=size(A); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if A(k,k)==0 return endm=A(i,k)/A(k,k); for j=k+1:nA(i,j)=A(i,j)-m*A(k,j); endb(i)=b(i)-m*b(k); enddet=det*A(k,k); %计算行列式enddet=det*A(n,n);for k=n:-1:1 %回代求解for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end2. 列主元Gauss消去法：function x=detGauss(A,b)% Gauss列主元消去法[n,m]=size(A);nb=length(b);det=1; %存储行列式值x=zeros(n,1);for k=1:n-1amax=0; %选主元for i=k:nif abs(A(i,k))>amaxamax=abs(A(i,k));r=i;endendif amax<1e-10return;endif r>k %交换两行for j=k:nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=-det;endfor i=k+1:n %进行消元m=A(i,k)/A(k,k);for j=k+1:nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A(n,n);for k=n:-1:1 %回代求解for j=k+1:nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);end矩阵A和b的构造clc;clear;n=10;%n=84;A=eye(n)*6+diag(ones(1,n-1)*8,-1)+diag(ones(1,n-1),1); b=[7,15*ones(1,n-2),14]';计算结果：(1)n=10时Gauss消元法>>x=DelGauss(A,b)x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000列主元Gauss消去法>>x=detGauss(A,b)x =1111111111(2) n=84时Gauss消元法>>x=DelGauss(A,b) x =1.0e+008 *0.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0002 -0.0003 0.0007 -0.0013 0.0026 -0.0052 0.0105 -0.0209 0.0419 -0.0836 0.16650.6501-1.25822.3487-4.02635.3684列主元Gauss消去法>>x=detGauss(A,b) x =1.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.00001.0000 1.0000 1.00001.00001.0000 1.0000结果分析由上述实验结果可知，对于n=10采用Gauss 消去法和Gauss 列主元消去法得到的实验结果是相同的，而对于n=84，Gauss 消去法所得到的结果是错误的，Gauss 列主元消去法得到的结果是正确的。
大连理工大学-矩阵大作业
（1）从大到小的顺序的计算程序：function y=snd(n)format longy=0;if n<2disp('请输入大于1的数！')elses=0;i=2;while i<=ns=single(s+(1/(i^2-1))); i=i+1;endy=s;end（2）从小到大的顺序的计算程序：function y=snx(n)format longy=0;if n<2disp('请输入大于1的数！')elses=0;i=n;while 1s=single(s+(1/(i^2-1)));i=i-1;if i==1breakendendy=s;end（3）按两种顺序分别计算并指出有效位数（编制程序时用单精度）S的计算结果：①210S的计算结果：②410S的计算结果：③610计算时的有效位数为七位数。
① 秦九昭算法计算程序：function y=qjz(a,x) j=3; i=size(a,2); switch i case 1 y=a(1); case 2y=a(1)*x+a(2); otherwise p=a(1)*x+a(2); while j<=i p=p*x+a(j); j=j+1; end y=p; end② 计算在点23的值。
计算结果如下：当23=x 时()86652=x f 。
①Gauss法计算程序和结果：程序：A(1,:)=[31,-13,0,0,0,-10,0,0,0]; A(2,:)=[-13,35,-9,0,-11,0,0,0,0]; A(3,:)=[0,-9,31,-10,0,0,0,0,0];A(4,:)=[0,0,-10,79,-30,0,0,0,-9]; A(5,:)=[0,0,0,-30,57,-7,0,-5,0]; A(6,:)=[0,0,0,0,-7,47,-30,0,0];A(7,:)=[0,0,0,0,0,-30,41,0,0];A(8,:)=[0,0,0,0,-5,0,0,27,-2];A(9,:)=[0,0,0,-9,0,0,0,-2,29];B=[-15;27;-23;0;-20;12;-7;7;10]; [a,b]=gauss(A,B);j=size(a,2);while j>=1k=j+1;s=b(j);while k<=9s=s-x(k)*a(j,k);k=k+1;endx(j)=s/a(j,j);j=j-1;enddisp(x)function [x,y]=gauss(a,b)num_i=size(a,1);j=1;while j<=(num_i-1)i=j+1;while i<=num_ir=a(i,j)/a(j,j);a(i,:)=a(i,:)-r*a(j,:);b(i,:)=b(i,:)-r*b(j,:);i=i+1;endj=j+1;endx=a;y=b;运行的结果为：()T 289.0-345.0.0=。
大连理工大学-2015年矩阵上机编程作业
} //Ax=b 的解答，A 为上三角矩阵 void SolveOne(double a[9][10],int m,int n){
int j=n-2; double answer[9]; for (int s=m-1;s>=0;s--){
answer[j]=a[s][n-1]/a[s][j]; for(int i=0;i<=s-1;i++){
} } }
} // 发现绝对值最大的元素所在一行和当前的一行做交换 void change(double a[9][10],int m,int n,int nowi,int nowj){
int record=nowi;//记录该和哪一行做交换 bool flag=true; //标志是否需要做交换 double max=a[nowi][nowj]; for(int i=nowi+1;i<m-1;i++){
}
//Gauss 消元解法 //先转化为上三角矩阵 void Gauss(double a[9][10],int m,int n){
double r[9]; //r 为每次的倍数 int k=0; for (int nowi=0,nowj=0;nowi<m-1 && nowj<n-1;nowi++,nowj++){
1. 设
, 其精确值为
.
(1) 编制按从大到小的顺序
, 计算 的通
用程序 (2) 编制按从小到大的顺序
, 计算
的通用程序
(3) 按两种顺序分别计算
并指出有效位数（编制程序时用单精度）
大连理工大学 矩阵与数值分析 第5章5.2.4Hermite与分段低次2017(春)
=
sin
⎛ ⎜⎝
π
2
⎞ ⎟⎠
=
1
,
f
′ ⎛⎜⎝
π
2
⎞ ⎟⎠
=
cos
⎛ ⎜⎝
π
2
⎞ ⎟⎠
=
0
,
H3 (x) =
f
⎛ ⎜⎝
π
4
⎞ ⎟⎠
⎛ ⎜ ⎜1 ⎜
−
2
x
π
−π
4
−π
⎝ 42
⎞ ⎟ ⎟ ⎟
⋅
⎛ ⎜ ⎜ ⎜
x
π
−π
2
−π
⎠⎝4 2
⎞2 ⎟ ⎟ ⎟ ⎠
+
f
′
⎛ ⎜⎝
π
4
⎞ ⎟⎠
⎛ ⎜⎝
x
−
π
4
⎞ ⎟⎠
⎛ ⎜ ⋅⎜ ⎜
二点三次Hermite插值多项
设已知函数 f(x) 在 2 个互异点 x1, x2 处的函数值和导数值：
f ( x1 ) , f ′( x1 ) ; f ( x2 ) , f ′( x2 ) ;
构造一个三次多项式 H3 ( x) = ax3 + bx2 + cx + d
使其满足插值条件
H3 ( x1 ) = f ( x1 ) , H3′ ( x1 ) = f ′( x1 ) ;
答案是否定的。原因在于插值节点增多导致插值多项式 的次数增高，而高次多项式的振荡次数增多有可能使插值多 项式在非节点处的误差变得很大。
例如，对于函数
f
(x
)
=
1
5 +x
2
在[-5,5]上构造等距节
xk
大连理工大学矩阵与数值分析上机作业13388
大连理工大学矩阵与数值分析上机作业课程名称：矩阵与数值分析研究生姓名：交作业日时间：2016 年12 月20日1.1程序：Clear all;n=input('请输入向量的长度n：')for i=1:n;v(i)=1/i;endY1=norm(v,1)Y2=norm(v,2)Y3=norm(v,inf)1.2结果n=10 Y1 =2.9290Y2 =1.2449Y3 =1n=100 Y1 =5.1874Y2 =1.2787Y3 =1n=1000 Y1 =7.4855Y2 =1.2822Y3 =1N=10000 Y1 =9.7876Y2 =1.2825Y3 =11.3 分析一范数逐渐递增，随着n的增加，范数的增加速度减小；二范数随着n的增加，逐渐趋于定值，无群范数都是1.2.1程序clear all;x(1)=-10^-15;dx=10^-18;L=2*10^3;for i=1:Ly1(i)=log(1+x(i))/x(i); d=1+x(i);if d == 1y2(i)=1;elsey2(i)=log(d)/(d-1);endx(i+1)=x(i)+dx;endx=x(1:length(x)-1);plot(x,y1,'r');hold onplot(x,y2);2.2 结果2.3 分析红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大小始终为1，这主要是由于大数加小数的原因。
第3题3.1 程序clear all;A=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512];x=1.95:0.005:2.05;for i=1:length(x);y1(i)=f(A,x(i));y2(i)=(x(i)-2)^9;endfigure(3);plot(x,y1);hold on;plot(x,y2,'r');F.m文件function y=f(A,x) y=A(1);for i=2:length(A); y=x*y+A(i); end;3.2 结果第4题4.1 程序clear all;n=input('请输入向量的长度n：')A=2*eye(n)-tril(ones(n,n),0);for i=1:nA(i,n)=1;endn=length(A);U=A;e=eye(n);for i=1:n-1[max_data,max_index]=max(abs(U(i:n,i))); e0=eye(n);max_index=max_index+i-1;U=e0*U;e1=eye(n);for j=i+1:ne1(j,i)=-U(j,i)/U(i,i);endU=e1*U;P{i}=e0;%把变换矩阵存到P中L{i}=e1;e=e1*e0*e;endfor k=1:n-2Ldot{k}=L{k};for i=k+1:n-1Ldot{k}=P{i}*Ldot{k}*P{i};endendLdot{n-1}=L{n-1};LL=eye(n);PP=eye(n);for i=1:n-1PP=P{i}*PP;LL=Ldot{i}*LL;endb=ones(n,2);b=e*b; %解方程x=zeros(n,1);x(n)=b(n)/U(n,n);for i=n-1:-1:1x(i)=(b(i)-U(i,:)*x)/U(i,i);endX=U^-1*e^-1*eye(n);%计算逆矩阵AN=X';result2{n-4,1}=AN;result1{n-4,1}=x;fprintf('%d:\n',n)fprintf('%d ',AN);4.2 结果n=51.0625 -0.875 -0.75 -0.5 -0.06250.0625 1.125 -0.75 -0.5 -0.06250.0625 0.125 1.25 -0.5 -0.06250.0625 0.125 0.25 1.5 -0.0625-0.0625 -0.125 -0.25 -0.5 0.0625n=101.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 1.125 -0.75 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 1.25 -0.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 0.0625 0.125 0.25 1.5 -0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625 -0.0625 -0.125 -0.25 -0.5 0.0625同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。
大连理工大学矩阵与数值分析MATLAB上机实验
二、解线性方程组 1．分别 Jacobi 迭代法和 Gauss-Seidel 迭代法求解线性方程组
3 1 0 0 x1 1 1 3 1 0 x2 0 , 0 1 2 1 x3 0 0 0 1 3 x4 0
Gauss 列主元消去法程序:
clc; clear; format long a=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5]; %系数矩阵 b=[12;6;23;16]; [n,m]=size(a); nb=length(b); det=1; for k=1:n-1 amax=0; for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k)); r=i; end end if amax<1e-10 return; end if r>k for j=k:n z=a(k,j); a(k,j)=a(r,j); a(r,j)=z; end z=b(k);
从小到大求和程序计算结果：
N 100 10000 1000000 从小到大求和程序得 到的 SN 0.497512437810945 0.499975001249937 0.499999750000134 真实值������������ = ������ 2������ + 1 0.497512437810945 0.499975001249937 0.499999750000125 计算值有效位数 15 15 13
8
2
1 dx x
复化梯形公式程序
clc; clear; format long syms t m=int(1/t,2,8); %真实值 a=2; b=8; n=300; h=(b-a)/n; sum=0; f=inline('1/x'); for i=1:n-1 sum=sum+f(a+i.*h); end T=h/2*(f(a)+2*sum+f(b))
大连理工大学矩阵与数值分析上机作业13388
共享知识分享快乐大连理工大学矩阵与数值分析上机作业课程名称：矩阵与数值分析研究生姓名：12 交作业日时间：日20 月年2016卑微如蝼蚁、坚强似大象．共享知识分享快乐第1题1.1程序：Clear ;all n=input('请输入向量的长度n：') for i=1:n;v(i)=1/i;endY1=norm(v,1)Y2=norm(v,2)Y3=norm(v,inf)1.2结果n=10 Y1 =2.9290Y2 =1.2449Y3 =1n=100 Y1 =5.1874Y2 =1.2787Y3 =1n=1000 Y1 =7.4855Y2 =1.2822Y3 =1N=10000 Y1 =9.7876Y2 =1.2825Y3 =11.3 分析一范数逐渐递增，随着n的增加，范数的增加速度减小；二范数随着n的增加，逐渐趋于定值，无群范数都是1.卑微如蝼蚁、坚强似大象．共享知识分享快乐第2题2.1程序;clear all x(1)=-10^-15;dx=10^-18;L=2*10^3; i=1:L fory1(i)=log(1+x(i))/x(i); d=1+x(i); d == 1ify2(i)=1;elsey2(i)=log(d)/(d-1);endx(i+1)=x(i)+dx;end x=x(1:length(x)-1););'r'plot(x,y1,on holdplot(x,y2);卑微如蝼蚁、坚强似大象．共享知识分享快乐2.2 结果2.3 分析红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大小始终为1，这主要是由于大数加小数的原因。
第3题3.1 程序;clear all A=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512];x=1.95:0.005:2.05; i=1:length(x);for y1(i)=f(A,x(i)); y2(i)=(x(i)-2)^9;end figure(3);plot(x,y1);;on hold卑微如蝼蚁、坚强似大象．共享知识分享快乐);'r'plot(x,y2,F.m文件y=f(A,x)function y=A(1); i=2:length(A);for y=x*y+A(i);;end3.2 结果第4题卑微如蝼蚁、坚强似大象．共享知识分享快乐4.1 程序;clear all n=input('请输入向量的长度n：')A=2*eye(n)-tril(ones(n,n),0); i=1:n for A(i,n)=1;end n=length(A);U=A; e=eye(n);for i=1:n-1[max_data,max_index]=max(abs(U(i:n,i))); e0=eye(n);max_index=max_index+i-1; U=e0*U; e1=eye(n); j=i+1:n fore1(j,i)=-U(j,i)/U(i,i);endU=e1*U;中把变换矩阵存到P P{i}=e0;% L{i}=e1; e=e1*e0*e;endk=1:n-2for Ldot{k}=L{k}; i=k+1:n-1forLdot{k}=P{i}*Ldot{k}*P{i};endend Ldot{n-1}=L{n-1};LL=eye(n);PP=eye(n); i=1:n-1for PP=P{i}*PP;LL=Ldot{i}*LL;endb=ones(n,2);解方程 %b=e*b;x=zeros(n,1);x(n)=b(n)/U(n,n); i=n-1:-1:1for卑微如蝼蚁、坚强似大象．共享知识分享快乐x(i)=(b(i)-U(i,:)*x)/U(i,i);end计算逆矩阵%X=U^-1*e^-1*eye(n);AN=X'; result2{n-4,1}=AN;result1{n-4,1}=x;,n)'%d:\n'fprintf(fprintf('%d ',AN);4.2 结果n=51.0625 -0.875 -0.75 -0.5 -0.0625-0.0625 0.0625 -0.75 1.125 -0.5-0.0625 0.125 0.0625 1.25 -0.5-0.0625 0.1250.25 0.06251.50.0625-0.5-0.25-0.0625 -0.125n=101.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 -0.0625 1.125 0.0625 -0.75 -0.5 -0.5 0.0625 1.125 -0.75 -0.0625 -0.0625 0.0625 0.125 1.25 1.25 -0.0625 -0.5 0.0625 0.125 -0.5-0.0625 0.250.250.0625 0.1251.5 1.5 -0.0625 0.1250.06250.0625 -0.0625 -0.125 -0.25 0.0625 -0.5 -0.0625 -0.125 -0.25 -0.5 -0.0625 -0.75 1.0625 -0.5 -0.0625 -0.875 -0.5 -0.75 1.0625 -0.875 -0.0625 -0.5 0.0625 1.125 -0.5 0.0625 1.125 -0.75 -0.0625 -0.75 1.25 0.125 0.0625 -0.0625 -0.0625 -0.5 -0.5 0.0625 0.125 1.250.25-0.0625 -0.0625 1.50.1250.0625 0.0625 0.250.1251.5-0.0625 -0.125 -0.25 0.0625-0.5 0.0625 -0.0625 -0.125 -0.25-0.5同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2013级工科硕士研究生《矩阵与数值分析》课程数值实验报告大连理工大学Dalian University of Technology一、设622101NNjSj==-∑，分别编制从小到大和从大到小的顺序程序分别计算100001000000,S S并指出两种方法计算结果的有效位数。
程序代码：从小到大：function f=s(N); %定义函数sf=0; %初始值为0for j=N:-1:3 %j从3到n循环（从小到大） ft=1000000/(j^2-1); %Sjf=f+ft; %SNend从大到小：function f=s(N); %定义函数sf=0; %初始值为0for j=N:-1:3 %j从3到n循环（从小到大） ft=1000000/(j^2-1); %Sjf=f+ft; %SNend执行结果：从小到大：s(10000)ans =4.16566671666167e+05s(1000000)ans =4.166656666671731e+05 有效数字：16,16 从大到小： s(10000) ans =4.165666716661668e+05s(1000000) ans =4.166656666671667e+05 有效数字：16,16 分析：小数和大数相加时，按照从大到小的顺序和按照从小到大的顺序得出的结果不同，前者由于舍入误差的影响而使结果不准确，所以应避免大数吃小数的现象。
二、解线性方程组1．分别利用Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组Ax b =，其中常向量为()21n -维随机生成的列向量，系数矩阵A 具有如下形式1111111122n n n n n n n n T I I I A I I T I --------+-⎛⎫ ⎪- ⎪= ⎪- ⎪-+⎝⎭， 其中1211112n T --⎛⎫⎪- ⎪= ⎪- ⎪-⎝⎭为1n -阶矩阵，1n I -为1n -阶单位矩阵，迭代法计算停止的条件为：101210k k x x -+-<，给出10,100,1000n =时的不同迭代步数． 程序代码：Jacobi迭代法function [u,l,b,x,kk,delta,A]=Ja(n)t=2.*eye(n-1);s=-eye(n-2);z=zeros(n-2,1); %生成n-1行1列零向量z1=zeros(1,n-1);s1=[z s];s1=[s1;z1];s2=[s z];s2=[z1;s2];T=t+s1+s2; %生成Ti=-eye((n-1)*(n-2));a1=[zeros(((n-1)*(n-2)),n-1) i];a2=[a1;zeros(n-1,((n-1)*(n-1)))];b1=[i zeros(((n-1)*(n-2)),n-1)];b2=[zeros(n-1,((n-1)*(n-1)));b1];d=4.*eye((n-1)*(n-1)); %生成Db=round(randn((n-1)*(n-1),1)); %生成bA0=[];for k=1:1:n-1x=T+2.*eye(n-1);A0=blkdiag(A0,x); %生成AendA=A0+a2+b2; %生成Au=-(triu(A)-d); %生成Ul=-(tril(A)-d); %生成Lx=randn((n-1)*(n-1),1); %生成初始x%%以下为计算迭代部分B=inv(d)*(l+u);ff=inv(d)*b;kk=1;derta=10;delta=[];x0=x;while derta>10^-10&(kk<800);x0=x;x=B*x+ff;derta=norm(x-x0);kk=kk+1;delta=[delta derta];endplot(delta);End执行结果：运行结果：[u,l,b,x,kk,delta,A]=Ja(10)情况图1 迭代曲线收敛步数465 收敛值9.899081611039407e-11[u,l,b,x,kk,delta,A]=Ja(100)情况下图2 迭代步数48988 收敛终值9.997187164447103e-11 当n=1000时，矩阵规模超出内存。
程序代码：Gauss-Seidel迭代法function kk=Gs(n)t=2.*eye(n-1);s=-eye(n-2);z=zeros(n-2,1); %生成n-1行1列零向量z1=zeros(1,n-1);s1=[z s];s1=[s1;z1];s2=[s z];s2=[z1;s2];T=t+s1+s2; %生成Ti=-eye((n-1)*(n-2));a1=[zeros(((n-1)*(n-2)),n-1) i];a2=[a1;zeros(n-1,((n-1)*(n-1)))];b1=[i zeros(((n-1)*(n-2)),n-1)];b2=[zeros(n-1,((n-1)*(n-1)));b1];d=4.*eye((n-1)*(n-1)); %生成Db=round(randn((n-1)*(n-1),1)); %生成bA0=[];for k=1:1:n-1x=T+2.*eye(n-1);A0=blkdiag(A0,x); %生成AendA=A0+a2+b2; %生成Au=-(triu(A)-d); %生成Ul=-(tril(A)-d); %生成Lx=randn((n-1)*(n-1),1); %生成初始x%%以下为Gauss迭代法部分B=(d-l)\u;F=(d-l)\b;x0=x;derta=10;kk=0;delta=[];while derta>10^-10&(kk<800);x0=x;x=B*x+F;derta=norm(x-x0);kk=kk+1;delta=[delta derta];endplot(delta);End执行结果： f=Gs(10)图3 G-S 法n=10 迭代步数220 收敛终值9.623878789391471e-11f=Gs(100) 分析：Gauss-Seidel 迭代法是Jacobi 迭代法的改进，求解该线性方程组，Gauss-Seidel 迭代法和Jacobi 迭代法都收敛，并且Gauss-Seidel 迭代法收敛速度比Jacobi 迭代法快。
但是具体的问题要具体分析。
2. 用Gauss 列主元消去法、QR 方法求解如下方程组：12341211025034.17921281218x x x x -⎛⎫⎛⎫⎛⎫ ⎪⎪ ⎪ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎪ ⎪---⎝⎭⎝⎭⎝⎭程序代码：Gauss列主元消去法function x=Gsl(A,b)%%输入系数矩阵A和常数列向量b，高斯列主元法Ab=[A b]; %构成增广矩阵[m,n]=size(Ab);for i=1:1:mg=tril(Ab);g=abs(g);mx=find(g(:,i)==max(g(:,i))); %找到第m列最大值所在行数为mxAb([i mx],:)=Ab([mx i],:); %将第m行与mx行交换for k=(i+1):1:mxb=-Ab(k,i)/Ab(i,i);xc=xb*Ab(i,:);Ab(k,:)=Ab(k,:)+xc;endend%%至此高斯列主元结束，下面为回带法x=zeros(m,1);%初始化解向量xx=x+1;A=Ab(:,1:m);b=Ab(:,n);for i=m:-1:1%回带法求x向量x(i,1)=0;format long; %取长精度x(i,1)=(b(i,1)-A(i,:)*x)./A(i,i);endEnd执行结果：x=（-1 0 1 2）T程序代码：QR 法function [Q,R,x,y]=qrf(A,b) [Q,R]=qr(A);%求A 的QR 分解矩阵 y=Q\b;%利用回带法求解 x=R\y; End 执行结果：x=（-1 0 1 2）Ty = (5.2590 -13.6412 -3.1085 0.7734)T分析：从上边的编程可以看出Gauss 列主元消去法与QR 方法求解的结果相同。
三、非线性方程的迭代解法1．求方程()222sin ln 160x f x e x x x =++--=的根，迭代停止的条件为：10110k k x x -+-<；程序代码：%% 二分法迭代求方程的根，迭代区间选为[1,3] clear all; syms x;f=exp(x)+2*x^2+2*sin(x)-log(x)-16; %函数f(x)的表达式 i=0; %二分次数记数 a=1; %求根区间左端 b=3; %求根区间右端 e=1e-10;fa=subs(f,x,a); fb=subs(f,x,b);c=(a+b)/2; %计算区间中点fc=subs(f,x,c);while abs(fc)>=e; %判断f(c)是否为零点if fa*fc>=0; %判断左侧区间是否有根fa=fc;a=c;elsefb=fc;b=c;endc=(a+b)/2;fc=subs(f,x,c);i=i+1;endfprintf('采用二分法迭代求解方程的根为：%.10f \n',c);fprintf('采用二分法迭代求解方程的根的迭代次数：%d \n',i); 执行结果：采用二分法迭代求解方程的根为：1.9629226832采用二分法迭代求解方程的根的迭代次数：352．利用Newton迭代法求多项式432--+-x x x x360311=的所有实零点，注意重根的问题。
程序代码：%% Newton迭代法求多项式的所有实根clear all;syms x;f=x^4-3*x^3-3*x^2+11*x-6; % 多项式f(x)fx=diff(f); %对fx求导fxx=diff(fx);%x0=-1.5; %选取初始点x0=0.5;%x0=2.5;f0=subs(f,x,x0); fx0=subs(fx,x,x0); e=1e-5; k=0;y=x0-f0/fx0;while abs(y-x0)>=e x0=y;f0=subs(f,x,x0); fx0=subs(fx,x,x0); y=x0-f0/fx0; k=k+1; endif (subs(fxx,x,y)~=0)&&(subs(fx,x,y)<=1e-4)fprintf('用Newton 迭代法所求得的多项式的根是2重根:\n'); fprintf('根为：%.4f \n,',y); fprintf('迭代次数为：%d \n',k); elsefprintf('用Newton 迭代法所求得的多项式的根是单根:\n'); fprintf('根为：%.4f \n,',y); fprintf('迭代次数为：%d \n',k); end执行结果：(1).用Newton 迭代法所求得的多项式的根是2重根:根为：1.0000 ,迭代次数为：15 (2). 用Newton 迭代法所求得的多项式的根是单根:根为：-2.0000 ,迭代次数为：5 (3). 用Newton 迭代法所求得的多项式的根是单根:根为：3.0000 ,迭代次数为：7分析：本程序利用牛顿法迭代法进行迭代计算，公式为)()('1k k k k x f x f mx x -=+, 另外牛顿迭代法的初始点应该在零点附近，否则可能迭代不收敛。