2009年大连理工大学矩阵与数值分析上机程序

大连理工线代上机实验

基本操作
四则运算、转置、求逆、求秩、求行列式、组合、化为行最简形、求特征值
常见任务
① 矩阵的赋值和其加、减、乘、除（求逆）命令； ② 矩阵化为最简行阶梯型的计算命令；[U0,ip]=rref(A) ③ 多元线性方程组MATLAB求解的几种方法；x=inv(A)*b, U=rref(A) ④ 行列式的几种计算机求解方法； D=det(A),[L,U]=lu(A);D=prod(diag(L)) ⑤ n个m维向量组的相关性及其秩的计算方法和命令; r=rank(A),U=rref(A) ⑥ 求欠定线性方程组的基础解系及超定方程解的MATLAB 命令；xb=null(A) ⑦ 矩阵的特征方程、特征根和特征向量的计算命令； f=poly(A);[P,D]=eig(A) ⑧ 化二次型为标准型的MATLAB命令；yTDy=xTAx; 其中 y=P-1x,
• • • •
-66.5556 25.6667 -18.7778 26.5556
• >>
例三、求秩
• • • • • • • • • • • • >> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >> r=rank(A); >> r= % = 计算机不显示r的值 ??? r= | Error: Expression or statement is incomplete or incorrect. >> rank(A) ans = 4 >> r r= 4 %不打；则计算机将显示rank(A)的值
例1 用直接解法求解下列线性方程组. 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; x=A\b

第一章矩阵与数值分析

以求解20阶线性方程组为例，如果用Cramer法则求解，在算法中的乘、除运算次数将大于
21！（约9.7×1020次）
使用每秒一亿次的串行计算机计算,完成运算耗时约30万年!
Cramer算法是“实际计算不了”的。为此，人们研究出著名的Gauss消去法，它的计算过程已作根本改进，使得上述例子的乘、除运算仅为3060次，这在任何一台电子计算机上都能很快完成。
这是由于
x a x a x a2
a
x
ax
a
(x a)2
a x
a
x
a
a
2
1
1 xa
x
a
2
a
a
是 x 的a平方项级，故可忽略不计。
a
相对误差也可正可负，其绝对值的上界叫做相对误差界
（限）。当绝对误差界为 e时a ，相对误差界取为
ea
相对误差界（限）
a
例1 已知 e 2.71828182 ,其近似值 a 2.718,求 a
1.2.1 误差来源与分类
用计算机解决科学计算问题时经常采用的处理方式是将连续的问题离散化、用有限代替无限等，并且用数值分析所处理的一些数据，不论是原始数据，还是最终结果，绝大多数都是近似的，因此在此过程中，误差无处不在．误差主要来源于以下四个方面：
实际问题
模型误差
计
算机科学计Βιβλιοθήκη 数学模型主要内容包括：
数值代数
Ax b
f x 0
数值逼近微分方程数值解法
f x f x
b
a
x
f
x
dx
u f t, u, ut0 u0
具体任务：一、构造在计算机上可行的有效算法。二、给出可靠的理论分析：进行误差分析，讨论数值算法的收敛性和数值稳定性。

矩阵与数值分析上机实验题及程序

1.给定n 阶方程组Ax b =，其中6186186186A ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ ，7151514b ⎛⎫ ⎪ ⎪ ⎪= ⎪⎪⎪⎝⎭则方程组有解(1,1,,1)Tx = 。

对10n =和84n =，分别用Gauss 消去法和列主元消去法解方程组，并比较计算结果。

Gauss 消去法：Matlab 编程（建立GS.m 文件）：function x=GS(n) A=[];b=[]; for i=1:n-1 A(i,i)=6; A(i,i+1)=1; A(i+1,i)=8; b(i)=15; endA(n,n)=6;b(1)=7;b(n)=14;b=b'; for k=1:n-1 for i=k+1:nm(i,k)=A(i,k)/A(k,k);A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); b(i)=b(i)-m(i,k)*b(k); end endb(n)=b(n)/A(n,n); for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i); end clear x; x=b;disp( 'AX=b 的解x 是') end计算结果：在matlab 命令框里输出GS （10）得： >> GS(10)AX=b 的解x 是 ans =1.0000 1.00001.00001.00001.00001.00001.00001.00001.0000在matlab命令框里输出GS（84）得：>> GS(84)AX=b的解x是ans =1.0e+008 *0.0000………0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00010.0002-0.00030.0007-0.00130.0026-0.00520.0105-0.02090.0419-0.08360.1665-0.3303-1.25822.3487-4.02635.3684列主元消去法：Matlab编程（建立GLZX.m文件）：function x=GLZX(n)A=[];b=[];eps=10^-2;for i=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8;b(i)=15;endA(n,n)=6;b(1)=7;b(n)=14;b=b';for k=1:n-1[mainElement,index]=max(abs(A(k:n,k)));index=index+k-1;%indexif abs(mainElement)<epsdisp('列元素太小！！');break;elseif index>ktemp=A(k,:);temp1=b(k);A(k,:)=A(index,:);b(k)=b(index);A(index,:)=temp;b(index)=temp1;endfor i=k+1:nm(i,k)=A(i,k)/A(k,k);%A(k,k) ;A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n);b(i)=b(i)-m(i,k)*b(k);endendb(n)=b(n)/A(n,n);for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n).*b(i+1:n)'))/A(i,i); endclear x;x=b;disp('AX=b的解x是')end在matlab命令框里输出GLZX（10）得：>> GLZX(10)AX=b的解x是ans =1111111111在matlab命令框里输出GLZX（84）得：>> GLZX(84)AX=b的解x是ans =1.00001.00001.00001.00001.00001.0000………1.00001.00001.00001.00001.00001.00001.00001.0000分析：n较小时，两种方法均能得到正确解，当n较大后，Gauss消去法计算结果严重偏离准确值成为错解，列主元消去法依然能得到正确解。

大连理工大学《矩阵与数值分析》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 。

大连理工大学-矩阵与数值分析第一章上

–对C的类型系统改进和扩充(更安全)
–支持面向对象
C++保持与C兼容(快速普及) C++不是纯粹的面向对象的语言

33
1.2 程序的编译过程
34
1.3 C++的词法记号
关键字各种常量操作符标识符分隔符

35
1.4 C++程序的结构
#include <iostream.h> int main() { cout<<”this is the start of something wonderful!”; cout<<endl; cout<<”And now we can say even more!”; return 0; }
Sub1
Sub2
….
Subn
各子流程实现----函数化 Func1 Func2 …. Funcn
根据系统的流程组建软件，通过函数的调用实现
17
面向对象思想
问题域 (Domain) 以问题域中的事物为中心思考问题 Object1 Object2
….
Objectn
对象归类----抽象化 Class1 Class2 …. Classn
返回类型
{ 函数体； }
50
函数名（形式参数1，形式参数2，。。。，形式参数
3.2 参数的传递
值调用
#include <iostream.h> double Volume(double radius,double height); int main() { double v; v=Volume(3.0,3.0); cout<<"Volume="<<v<<endl; return 0; } double Volume(double radius,double height) { double result=3.14 * radius * radius * height; return result; }

大连理工大学矩阵与数值分析上机作业代码

则方程组有解 x = (1,1,⋯ ,1) 。对 n = 10 和 n = 84 ，分别用 Gauss 消去法和列主元消去法解
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;结果三次函数的均方误差最小，拟合的最好。

大连理工大学矩阵与数值分析上机作业

s=s+abs(x(i));
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=zeros(n,1); %为列向量 x 预分配存储空间 y=1:n; %定义行向量 y y=y'; %把行向量 y 改成列向量 for i=1:n
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 为输入参数

[理学]矩阵与数值分析-第1章li-Chapter1

什么是有效算法？
考察，线性方程组的解法
⎧a11 x1 + a12 x2 + L + a1n xn = b1 ⎪a x + a x + L + a x = b ⎪ 21 1 22 2 2n n 2 ⎨ M ⎪ ⎪ ⎩an1 x1 + an 2 x2 + L + ann xn = bn
早在18世纪Cramer已给出了求解法则：
2001. 数值线性代数，徐树方、高立、张平文，北京大学出版社，2000. 数值逼近，王仁宏，高等教育出版社，1999. 数值逼近方法，南京大学数学系计算数学专业编，科学出版社，1978. 微分方程数值解法，李荣华、冯果忱，高等教育出版社， 1996. 微分方程数值方法，胡健伟、汤怀民，科学出版社，1999. 矩阵分析引论，罗家洪、方卫东，华南理工大学出版社， 2006. 矩阵分析，同济大学应用数学系，同济大学出版社，2005.
本课程的成绩考核标准
1、平时的课后作业 2、数值试验报告 (Matlab,C) 3、期末考试 ≈ 70%
课程网站 /numerical/
≈ 30%
第1章
绪
论
1.1 计算机科学计算研究对象与特点
科学计算、理论计算和实验并列为三大科学方法。现代意义下的计算数学主要研究在计算机上计算的有效算法及其相关理论，从而使它成为一门新学科——科学计算。
矩阵与数值分析
大连理工大学工科硕士研究生基础学位课程
李崇君
(主讲) 作者：张宏伟、金光日、李崇君大连理工大学数学科学学院
从一个小例子开始
问题: 在一个正方形的桌面上(边长为a), 分别在四个角有四只小虫,它们同时向着逆时针方向的另一只小虫移动, 求它们的运动轨迹.

大工矩阵上机题答案

上机题一．（1）s=0;for j=2:100;s=s+1/(j^2-1); endss =0.7400s=0;for j=100:-1:2;s=s+1/(j^2-1); endss =0.7400（2）s=0;for j=2：10000;s=s+1/(j^2-1); endss =0.7499for j=10000:-1:2;s=s+1/(j^2-1); endss =0.7499（3）s=0;for j=2：1000000;s=s+1/(j^2-1); endss =0.7500s=0;for j=1000000:-1:2; s=s+1/(j^2-1); endss =0.7500二1、Jacobi 迭代法算法：对于线性方程组Ax=b ，如果A 为非奇异方程，则可将A 分解为：A=D-L-U 其中D 为对角阵，其元素为A 的对角元素，L 与U 为A 的下三角阵和上三角阵。

于是Ax=b 化为：111()k k x D L U x D b --+=++，其中1()J B D L U -=+，1f D b -=。

程序：function x=jacobi(A,b,x0)A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 ;0; 0; 0];x0=[0;0 ;0 ;0];D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=D\(L+U);f=D\b;x=B*x0+f;n=1;while norm(x-x0,2)>=1.0e-6x0=x;x=B*x0+f;n=n+1;endfprintf('迭代次数为：')nfprintf('方程组的解为：')计算结果：迭代次数为：n =60方程组的解为：ans =0.80000.60000.40000.2000.Gauss-seidel 迭代法function x=Gaussseidel(A,b,x0)A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 ;0; 0; 0];x0=[0;0 ;0 ;0];D=diag(diag(A));U=-triu(A,1);L=-tril(A,-1);B=(D-L)\U;f=(D-L)\b;x=B*x0+f;n=1;while norm(x-x0,2)>=1.0e-6x0=x;x=B*x0+f;n=n+1;endfprintf('迭代次数')nfprintf('方程组的解为')迭代次数为：n =31方程组的解为：ans =0.80000.60000.40000.20002.用Gauss列主元消去法算法：Gauss列主元消去法是在Gauss消去法中增加选主元的过程，即在第k步（k=1,2,3,…）消元时，首先在第k列主对角元以下（含对角元）元素中挑选绝对值最大的数（即为列主元），并通过初等行变换，使得该数位于主对角线上，然后再继续消元。

大连理工大学矩阵与数值分析上机作业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时的结果，这里就不罗列了。

大连理工大学程名松矩阵与数值分析计算方法-第3章

1 (1 − ρ ( A ) ) > 0 一定存在 2 一种相容的矩阵范数 ⋅ ，使得 A ≤ ρ ( A) + ε 。
充分性根据定理2.8，对于 ε =
又根据相容矩阵范数的性质，再注意到上述关系式有 1 ρ ( A ) + ε = (1 + ρ ( A ) ) < 1 2 那么
A k ≤ A ≤ ( ρ ( A) + ε ) ≤ q k < 1
k→∞
k→∞
1⎞ lim ⎛ 1 + ⎟ k →∞ ⎜ ⎝ k⎠
k
lim sin k k →∞ k lim
lim
k →∞
k →∞
e− k
k
k
⎞ ⎛e 0 ⎞ ⎟ ⎜1 0 ⎟= A ⎟ =⎜ ⎟ ⎟ ⎜1 1 ⎟ ⎠ ⎠ ⎝
由矩阵序列极限的定义可以看出，矩阵序列收敛的性质和数列收敛性质相似。由定义可见，C m × n 中的矩阵序列的收敛相当于mn个数列同时收敛。因此可以用初等分析的方法来研究它。
∞ ∞
并且则
lim Ak = A ， k →∞
lim Bk = B
k →∞
lim A k B k = AB
k →∞
证由
Ak Bk − AB = Ak Bk − Ak B + Ak B − AB
≤ B ⋅ Ak − A + Ak ⋅ Bk − B
由定理1和推论可知，结论成立。
性质3
n n lim Ak = A 并且设 {Ak }k=1∈C × 中的矩阵序列，
∞
⎛ 0.1 0.7⎞ 。由于 A ∞ = 0.9 < 1 ，故计算 ∑ A ，其中 A = ⎜ ⎟ k =0 ⎝ 0.3 0.6⎠ ∞ k ρ ( A ) < 1，从而 ∑ A 收敛，且有