当前位置：文档之家 › 大连理工大学矩阵与数值分析上机作业

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

  • 格式：docx
  • 大小：449.86 KB
  • 文档页数：34

下载文档原格式

  / 34

合集下载

大连理工大学矩阵分析网上作业

大连理工大学矩阵分析网上作业

2010级工科硕士研究生《矩阵与数值分析》课程数值实验题目 一、设2211N N j S j ==−∑，分别编制从小到大和从大到小的顺序程序计算 100100001000000,,S S S ，并指出有效位数。

二、解线性方程组1．分别Jacobi 迭代法和Gauss ‐Seidel 迭代法求解线性方程组12342100112100,0121000120x x x x −−⎛⎞⎛⎞⎛⎞⎜⎟⎜⎟⎜⎟−⎜⎟⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟−⎜⎟⎜⎟⎜⎟−⎝⎠⎝⎠⎝⎠ 迭代法计算停止的条件为：6)()1(3110max −+≤≤<−k j k j j x x ．2. 用Gauss 列主元消去法、QR 方法求解如下方程组：12341212425327.2235113230x x x x ⎛⎞⎛⎞⎛⎞⎜⎟⎜⎟⎜⎟−⎜⎟⎜⎟⎜⎟=⎜⎟⎜⎟⎜⎟−−−⎜⎟⎜⎟⎜⎟⎝⎠⎝⎠⎝⎠ 三、非线性方程的迭代解法1．用Newton 迭代法求方程()22cos 60x x f x e x −=++−= 的根，计算停止的条件为：6110−+<−k k x x ；2．利用Newton 迭代法求多项式 43210.565.48 2.795.954=10x x x x −+−+的所有实零点，注意重根的问题。

四、数值积分分别用复化梯形公式和Romberg 公式计算积分821dx x∫要求误差不超过510−，并给出节点个数。

五、插值与逼近1．给定[]1,1−上的函数()22511xx f +=，请做如下的插值逼近： ⑴ 构造等距节点分别取5=n ，8=n ，10=n 的Lagrange 插值多项式；⑵ 构造分段线性取10=n 的Lagrange 插值多项式；⑶取Chebyshev 多项式()()x n x T n arccos cos ⋅=的零点： πnk x k 212cos−=，n k ,,1"= 作插值节点构造10=n 的插值多项式 ()x f 和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。

大连理工大学《矩阵与数值分析》学习指导与课后参考答案第三章、逐次逼近法

大连理工大学《矩阵与数值分析》学习指导与课后参考答案第三章、逐次逼近法

第三章 逐次逼近法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 法都收敛。

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

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

–对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;结果三次函数的均方误差最小，拟合的最好。

《矩阵与数值分析》上机大作业matlab

《矩阵与数值分析》上机大作业matlab

《矩阵与数值分析》上机大作业1.给定n 阶方程组Ax b =，其中6186186186A ⎛⎫ ⎪ ⎪ ⎪= ⎪ ⎪ ⎪⎝⎭ ，7151514b ⎛⎫ ⎪ ⎪ ⎪= ⎪⎪⎪⎝⎭则方程组有解(1,1,,1)Tx = 。

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

%产生三对角矩阵 n=84; %或n=10;A=zeros(n); b=zeros(1,n); for i=1:n-1A(i,i)=6;A(i,i+1)=1;A(i+1,i)=8; endA(n,n)=6;for i=2:n-1 b(1)=6; b(i)=15; b(n)=14; end Ab=[A b'];%Gauss 消元法for j=1:n-1 %按列循环 for k=j+1:n %消元Ab(k,:)=Ab(k,:)-Ab(j,:)*(Ab(k,j)/Ab(j,j)); end endx(n)=Ab(n,n+1)/Ab(n,n); for i=n-1:-1:1 %回代法求x for j=n:-1:i+1Ab(i,n+1)=Ab(i,n+1)-Ab(i,j)*x(j); endx(i)=Ab(i,n+1)/Ab(i,i); end（1）当n=10时，Gauss 消去法 Gauss 列主元消去法 x=1.000000000000000 x=1.000000000000000 1.000000000000000 1.0000000000000001.000000000000000 1.000000000000000 1.000000000000001 1.000000000000000 0.999999999999998 1.000000000000000 1.000000000000004 1.000000000000000 0.999999999999993 1.000000000000000 1.000000000000012 1.000000000000000 0.999999999999979 1.000000000000000 1.000000000000028 1.000000000000000（2） 当n=84时，Gauss 消去法的解是错解Columns 34 through 392147483649.00000 -4294967295.00000 8589934592.99999 -17179869182.9999 34359738368.9998Gauss 列主元消去法x 与x=(1,1…1)T 偏差不大 Columns 34 through 391.000000172108412 0.999999661246936 1.000000655651093 0.999998776117961 1.000002098083496综上，高斯列主元消去法可以避免小数作除数带来的误差，获得满意的数值解。

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

大连理工大学矩阵与数值分析上机作业
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上机作业

大连理工大学矩阵分析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 为输入参数

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

大连理工大学矩阵与数值分析上机作业.docx
functiony=QinJS(a,x)
务丫输出函数值
%a多项式系数，由高次到零次
显给定点、
n=length(a);
s=a(l);
fori=2:n
s=s*x+a (i);
end
Y=s；
计算
clear all;
clc;
・6 : 0・2:2・4 ;%>:=2的邻域
dispCx=2的邻域：f);x
a=[l -18 144 -672 2016 -4032 5376 -4608 2304 -512];
0.3333
0.2500
0.2000
0.1667
当n=7
方程精确解：
X =
1.0000
0.5000
0.3333
0.2500
0
0.1667
0.1429
利用
1.0000
0.5000
0.3333
0.2500
0.2000
0.1667
0.1429
利用
xPLU=
1.0000
0.5000
0.3333
0.2500
end
end
if
forj =1:m
c=U(izj)；
U(izj)=U(azj);
U(a, j)=c;
end
forj =1:m
c=P(izj)；
P(i,j)=P(a,j);
P(a,j)=c;
end
c=t (a)；
t (a)=t(i)；
t (i)=c；
end
U (i, i) =t (i)；
forj=i+l:m
z=z+L(izk)*U(kzj);

大连理工_2012矩阵与数值分析大作业

大连理工_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）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列主对角元以下（含对角元）元素中挑选绝对值最大的数（即为列主元），并通过初等行变换，使得该数位于主对角线上，然后再继续消元。

大连理工大学-2015年矩阵上机编程作业

大连理工大学-2015年矩阵上机编程作业
cout<<endl; } } }
} //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) 按两种顺序分别计算
并指出有效位数（编制程序时用单精度）

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

大连理工大学矩阵与数值分析上机作业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上机实验

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

大连理工大学矩阵与数值分析上机作业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章

大连理工大学程名松矩阵与数值分析计算方法-第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 收敛，且有
  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

矩阵与数值分析上机作业学校：大连理工大学学院：班级：姓名：学号：授课老师：注：编程语言Matlab程序：Norm.m函数function s=Norm(x,m)%求向量x的范数%m取1,2,inf分别表示1,2，无穷范数n=length(x);s=0;switch mcase 1 %1-范数for i=1:ns=s+abs(x(i));endcase 2 %2-范数for i=1:ns=s+x(i)^2;ends=sqrt(s);case inf %无穷-范数s=max(abs(x));end计算向量x，y的范数Test1.mclear all;clc;n1=10;n2=100;n3=1000;x1=1./[1:n1]';x2=1./[1:n2]';x3=1./[1:n3]'; y1=[1:n1]';y2=[1:n2]';y3=[1:n3]';disp('n=10时');disp('x的1-范数:');disp(Norm(x1,1));disp('x的2-范数:');disp(Norm(x1,2));disp('x的无穷-范数:');disp(Norm(x1,inf)); disp('y的1-范数:');disp(Norm(y1,1));disp('y的2-范数:');disp(Norm(y1,2));disp('y的无穷-范数:');disp(Norm(y1,inf)); disp('n=100时');disp('x的1-范数:');disp(Norm(x2,1));disp('x的2-范数:');disp(Norm(x2,2));disp('x的无穷-范数:');disp(Norm(x2,inf));disp('y的1-范数:');disp(Norm(y2,1));disp('y的2-范数:');disp(Norm(y2,2));disp('y的无穷-范数:');disp(Norm(y2,inf));disp('n=1000时');disp('x的1-范数:');disp(Norm(x3,1));disp('x的2-范数:');disp(Norm(x3,2));disp('x的无穷-范数:');disp(Norm(x3,inf));disp('y的1-范数:');disp(Norm(y3,1));disp('y的2-范数:');disp(Norm(y3,2));disp('y的无穷-范数:');disp(Norm(y3,inf));运行结果：n=10时x的1-范数:2.9290；x的2-范数:1.2449； x的无穷-范数:1y的1-范数:55； y的2-范数:19.6214； y的无穷-范数:10n=100时x的1-范数:5.1874；x的2-范数: 1.2787； x的无穷-范数:1y的1-范数:5050； y的2-范数:581.6786； y的无穷-范数:100n=1000时x的1-范数:7.4855； x的2-范数:1.2822； x的无穷-范数:1y的1-范数: 500500； y的2-范数:1.8271e+004；y的无穷-范数:1000程序Test2.mclear all;clc;n=100;%区间h=2*10^(-15)/n;%步长x=-10^(-15):h:10^(-15);%第一种原函数f1=zeros(1,n+1);for k=1:n+1if x(k)~=0f1(k)=log(1+x(k))/x(k);elsef1(k)=1;endendsubplot(2,1,1);plot(x,f1,'-r');axis([-10^(-15),10^(-15),-1,2]); legend('原图');%第二种算法f2=zeros(1,n+1);for k=1:n+1d=1+x(k);if(d~=1)f2(k)=log(d)/(d-1);elsef2(k)=1;endendsubplot(2,1,2);plot(x,f2,'-r');axis([-10^(-15),10^(-15),-1,2]); legend('第二种算法');运行结果：显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当]10,10[1515--∈x 时，x d +=1，计算机进行舍入造成d 恒等于1，结果函数值恒为1。

程序：秦九韶算法:QinJS.mfunction y=QinJS(a,x) %y 输出函数值%a 多项式系数，由高次到零次 %x 给定点n=length(a); s=a(1); for i=2:n s=s*x+a(i); end y=s;计算p(x):test3.mclear all ; clc;x=1.6:0.2:2.4;%x=2的邻域disp('x=2的邻域：');xa=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512]; p=zeros(1,5); for i=1:5p(i)=QinJS(a,x(i)); enddisp('相应多项式p 值：');p xk=1.95:0.01:20.5; nk=length(xk); pk=zeros(1,nk); k=1; for k=1:nkpk(k)=QinJS(a,xk(k)); endplot(xk,pk,'-r');xlabel('x');ylabel('p(x)');运行结果：x=2的邻域：x =1.6000 1.80002.0000 2.2000 2.4000相应多项式p值：p =1.0e-003 *-0.2621 -0.0005 0 0.0005 0.2621 p(x)在 x[1.95,20.5]上的图像程序：LU分解，LUDecom.mfunction [L,U]=LUDecom(A)%不带列主元的LU分解N = size(A);n = N(1);L=eye(n);U=zeros(n);for i=1:nU(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1); endfor i=2:nfor j=i:nz=0;for k=1:i-1z=z+L(i,k)*U(k,j);endU(i,j)=A(i,j)-z;endfor j=i+1:nz=0;for k=1:i-1z=z+L(j,k)*U(k,i);endL(j,i)=(A(j,i)-z)/U(i,i);endendPLU分解，PLUDecom.mfunction [P,L,U] =PLUDecom(A)%带列主元的LU分解[m,m]=size(A);U=A;P=eye(m);L=eye(m);for i=1:mfor j=i:mt(j)=U(j,i);for k=1:i-1t(j)=t(j)-U(j,k)*U(k,i);endenda=i;b=abs(t(i));for j=i+1:mif b<abs(t(j))b=abs(t(j));a=j;endendif a~=ifor j=1:mc=U(i,j);U(i,j)=U(a,j);U(a,j)=c;endfor j=1:mc=P(i,j);P(i,j)=P(a,j);P(a,j)=c;endc=t(a);t(a)=t(i);t(i)=c;endU(i,i)=t(i);for j=i+1:mU(j,i)=t(j)/t(i);endfor j=i+1:mfor k=1:i-1U(i,j)=U(i,j)-U(i,k)*U(k,j);endendendL=tril(U,-1)+eye(m);U=triu(U,0);（1）（2）程序：Test4.mclear all;clc;for n=5:30x=zeros(n,1);A=-ones(n);A(:,n)=ones(n,1);for i=1:nA(i,i)=1;for j=(i+1):(n-1)A(i,j)=0;endx(i)=1/i;enddisp('当n=');disp(n);disp('方程精确解:');xb=A*x; %系数bdisp('利用LU分解方程组的解:');[L,U]=LUDecom(A); %LU分解xLU=U\(L\b)disp('利用PLU分解方程组的解:');[P,L,U] =PLUDecom(A); %PLU分解xPLU=U\(L\(P\b))%求解A的逆矩阵disp('A的准确逆矩阵:');InvA=inv(A)InvAL=zeros(n); %利用LU分解求A的逆矩阵 I=eye(n);for i=1:nInvAL(:,i)=U\(L\I(:,i));enddisp('利用LU分解的A的逆矩阵:');InvALEnd运行结果：（1）只列出n=5,6,7的结果当n= 5方程精确解:x =1.00000.50000.33330.25000.2000利用LU分解方程组的解:xLU =1.00000.50000.33330.25000.2000利用PLU分解方程组的解:xPLU =1.00000.50000.33330.25000.2000当n=6方程精确解:x =1.00000.50000.33330.25000.20000.1667利用LU分解方程组的解: xLU =1.00000.50000.33330.25000.20000.1667利用PLU分解方程组的解: xPLU =1.00000.50000.33330.25000.20000.1667当n= 7方程精确解:x =1.00000.50000.33330.25000.20000.16670.1429利用LU分解方程组的解: xLU =1.00000.50000.33330.25000.20000.16670.1429利用PLU分解方程组的解:xPLU =1.00000.50000.33330.25000.20000.16670.1429（2）只列出n=5,6,7时A的逆矩阵的结果当n= 5A的准确逆矩阵:InvA =0.5000 -0.2500 -0.1250 -0.0625 -0.06250 0.5000 -0.2500 -0.1250 -0.12500 0 0.5000 -0.2500 -0.25000 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0625利用LU分解的A的逆矩阵:InvAL =0.5000 -0.2500 -0.1250 -0.0625 -0.06250 0.5000 -0.2500 -0.1250 -0.12500 0 0.5000 -0.2500 -0.25000 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0625当n= 6A的准确逆矩阵:InvA =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0625 0 0 0.5000 -0.2500 -0.1250 -0.1250 0 0 0 0.5000 -0.2500 -0.2500 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0313 利用LU分解的A的逆矩阵:InvAL =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.03130 0.5000 -0.2500 -0.1250 -0.0625 -0.06250 0 0.5000 -0.2500 -0.1250 -0.12500 0 0 0.5000 -0.2500 -0.25000 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0313当n= 7A的准确逆矩阵:InvA =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 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 0 0 0 0 0 0.5000 -0.5000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156利用LU分解的A的逆矩阵:InvAL =0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0156 -0.0156 0 0.5000 -0.2500 -0.1250 -0.0625 -0.0313 -0.0313 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 0 0 0 0 0 0.5000 -0.50000.5000 0.2500 0.1250 0.0625 0.0313 0.0156 0.0156程序：Cholesky分解：Cholesky.mfunction L=Cholesky(A)N = size(A);n = N(1);L=zeros(n);L(1,1)=sqrt(A(1,1));for i=2:nL(i,1)=A(i,1)/L(1,1);endfor j=2:ns1=0;for k=1:j-1s1=s1+L(j,k)^2;endL(j,j)=sqrt(A(j,j)-s1);for i=j+1:ns2=0;for k=1:j-1s2=s2+L(i,k)*L(j,k);endL(i,j)=(A(i,j)-s2)/L(j,j);endend计算Ax=b；Test5.mclear all;clc;for n=10:20A=zeros(n,n);b=zeros(n,1);for i=1:nfor j=1:nA(i,j)=1/(i+j-1);endb(i,1)=i;enddisp('n=');disp(n);disp('方程组原始解');x0=A\bdisp('利用Cholesky分解的方程组的解'); L=Cholesky(A)x=L'\(L\b)end运行结果：只列出了n=10,11的结果n=10方程组原始解x0 =1.0e+008 *-0.00000.0010-0.02330.2330-1.21083.5947-6.32336.5114-3.62330.8407利用Cholesky分解的方程组的解x =1.0e+008 *-0.00000.0010-0.02330.2330-1.21053.5939-6.32196.5100-3.62250.8405n=11方程组原始解x0 =1.0e+009 *0.0000-0.00020.0046-0.05670.3687-1.40393.2863-4.78694.2260-2.06850.4305利用Cholesky分解的方程组的解x =1.0e+009 *0.0000-0.00020.0046-0.05630.3668-1.39723.2716-4.76694.2094-2.06080.4290程序：（1）House.mfunction u=House(x)n=length(x);e1=eye(n,1);w=x-norm(x,2)*e1;u=w/norm(w,2);(2)Hou_A.mfunction HA=Hou_A(A)a1=A(:,1);n=length(a1);e1=eye(n,1);w=a1-norm(a1,2)*e1;u=w/norm(w,2);H=eye(n)-2*u*u'HA=H*A;(3)test6.mclear all;clc;A=[1 2 3 4;-1 2 sqrt(2) sqrt(3);-2 2 exp(1) pi;-sqrt(10) 2 -3 7;0 2 7 5/2];HA=Hou_A(A)运行结果：H =0.2500 -0.2500 -0.5000 -0.7906 0 -0.2500 0.9167 -0.1667 -0.2635 0 -0.5000 -0.1667 0.6667 -0.5270 0 -0.7906 -0.2635 -0.5270 0.1667 00 0 0 0 1.0000 HA =4.0000 -2.5811 1.4090 -6.53780.0000 0.4730 0.8839 -1.78050.0000 -1.0541 1.6576 -3.88360.0000 -2.8289 -4.6770 -4.10780 2.0000 7.0000 2.5000程序：Jacobi迭代：Jaccobi.mfunction [x,n]=Jaccobi(A,b,x0)%--·方程组系数阵A%--·方程组右端顶b%-- 初始值x0%--求解要求精确度eps%--迭代步数控制M%--·返回求得的解x%--·返回迭代步数nM=1000;eps=1.0e-5;D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵J=D\(L+U);f=D\b;x=J*x0+fn=1; %迭代次数err=norm(x-x0,inf)while(err>=eps)x0=x;x=J*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M)disp('Warning: 迭代次数太多，可能不收敛?');return;endendGauss_Seidel迭代：Gauss_Seidel.mfunction [x,n]=Gauss_Seidel(A,b,x0)%--Gauss-Seidel迭代法解线性方程组%--方程组系数阵 A%--方程组右端项 b%--初始值 x0%--求解要求的精确度 eps%--迭代步数控制 M%--返回求得的解 x%--返回迭代步数 neps=1.0e-5;M=10000;D=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+fn=1; %迭代次数err=norm(x-x0,inf)while(err>=eps)x0=x;x=G*x0+fn=n+1;err=norm(x-x0,inf)if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endend解方程组，test7.mclear all;clc;A=[5 -1 -3;-1 2 4;-3 4 15];b=[-2;1;10];disp('精确解');x=A\bdisp('迭代初始值');x0=[0;0;0]disp('Jacobi迭代过程：');[xj,nj]=Jaccobi(A,b,x0);disp('Jacobi最终迭代结果：');xjdisp('迭代次数');njdisp('Gauss-Seidel迭代过程：'); [xg,ng]=Gauss_Seidel(A,b,x0);disp('Gauss-Seidel最终迭代结果：'); xgdisp('迭代次数');ng运行结果：精确解x =-0.0820-1.80331.1311迭代初始值x0 =Jacobi迭代过程：x =-0.40000.50000.6667err =0.6667x =0.1000-1.03330.4533err =1.5333...x =-0.0820-1.80331.1311err =9.6603e-006Jacobi最终迭代结果：xj =-0.0820-1.80331.1311迭代次数nj =281Gauss-Seidel迭代过程：x =-0.40000.30000.5067err =0.5067x =-0.0360-0.53130.8012err =0.8313x =-0.0256-1.11510.9589err =0.5838...x =-0.0820-1.80331.1311err =9.4021e-006Gauss-Seidel最终迭代结果：xg =-0.0820-1.80331.1311迭代次数ng =20程序：Newton迭代法：Newtoniter.mfunction [x,iter,fvalue]=Newtoniter(f,df,x0,eps,maxiter) %牛顿法 x得到的近似解%iter迭代次数%fvalue函数在x处的值%f，df被求的非线性方程及导函数%x0初始值%eps 允许误差限%maxiter 最大迭代次数fvalue=subs(f,x0);dfvalue=subs(df,x0);for iter=1:maxiterx=x0-fvalue/dfvalueerr=abs(x-x0)x0=x;fvalue=subs(f,x0)dfvalue=subs(df,x0);if(err<eps)||(fvalue==0),break,endend弦截法：secant.mfunction [x,iter,fvalue]=secant(f,x0,x1,eps,maxiter)%弦截法 x得到的近似解%iter迭代次数%fvalue函数在x处的值%f被求的非线性方程%x0，x1初始值%eps 允许误差限%maxiter 最大迭代次数fvalue0=subs(f,x0);fvalue=subs(f,x1);for iter=1:maxiterx=x1-fvalue*(x1-x0)/(fvalue-fvalue0)err=abs(x-x1)x0=x1;x1=x;fvalue0=subs(f,x0);fvalue=subs(f,x1)if(err<eps)||(fvalue==0),break,endend求方程的实根：test8.mclear all；clc;syms xf=x.^3+2*x.^2+10*x-100;df=diff(f,x,1);eps=10e-6;maxiter=100;disp('Newton迭代初始值');xn1_0=0disp('Newton迭代结果');[xn1,iter_n1,fxn1]=Newtoniter(f,df,xn1_0,eps,maxiter) disp('Newton迭代初始值');xn2_0=5disp('Newton迭代结果');[xn2,iter_n2,fxn2]=Newtoniter(f,df,xn2_0,eps,maxiter) disp('弦截法初始值');xk1_0=0xk1_1=1disp('弦截法迭代结果');[xk1,iter_k1,fxk1]=secant(f,xk1_0,xk1_1,eps,maxiter) disp('弦截法初始值');xk2_0=5xk2_1=6disp('弦截法迭代结果');[xk2,iter_k2,fxk2]=secant(f,xk2_0,xk2_1,eps,maxiter)运行结果：Newton法结果:取两个不同初值0，5程序：二分法：resecm.mfunction [x,iter]=resecm(f,a,b,eps) %二分法 x 近似解%iter 迭代次数%f 求解的方程%a，b 求解区间%eps 允许误差限fa=subs(f,a);fb=subs(f,b);iter=0;if(fa==0)x=a;returnendif(fb==0)x=b;returnendwhile(abs(a-b)>=eps)mf=subs(f,(a+b)/2);if(mf==0)x=mf;n=n+1;returnendif(mf*fa<0)b=(a+b)/2;elsea=(a+b)/2;enditer=iter+1;endx=(a+b)/2;iter=iter+1;求方程的实根：test9.mclear all;clc;syms xf=exp(x).*cos(x)+2;a=0;a1=pi;a2=2*pi;a3=3*pi;b=4*pi;eps=10e-6;[x1,iter1]=resecm(f,a,a1,eps)[x2,iter2]=resecm(f,a1,a2,eps)[x3,iter3]=resecm(f,a2,a3,eps)[x4,iter4]=resecm(f,a3,b,eps)运行结果：[0,pi]区间的根x1 =1.8807；迭代次数iter1 = 20[pi,2*pi]区间的根x2 =4.6941；迭代次数iter2 =20[2*pi,3*pi]区间的根x3 =7.8548；迭代次数iter3 =20[3*pi,4*pi]区间的根x4 =10.9955；迭代次数iter4 =20程序：Newton插值：Newtominter.mfunction f=Newtominter(x,y,x0)%牛顿插值 x插值节点%y为对应的函数值%函数返回Newton插值多项式在x_0点的值fsyms t;if(length(x) == length(y))n = length(x);c(1:n) = 0.0;elsedisp('x和y的维数不相等！');return;endf = y(1);y1 = 0;l = 1;for(i=1:n-1)for(j=i+1:n)y1(j) = (y(j)-y(i))/(x(j)-x(i));endc(i) = y1(i+1);l = l*(t-x(i));f = f + c(i)*l;simplify(f);y = y1;if(i==n-1)if(nargin == 3) %如果3个参数则给出插值点的插值结果 f = subs(f,'t',x0);else%如果2个参数则直接给出插值多项式f = collect(f); %将插值多项式展开f = vpa(f, 6);endendend用等距节点做f(x)的Newton插值：test10.mn1=5;n2=10;n3=15;x0=[0:0.01:1];y0=sin(pi.*x0);x1=linspace(0,1,n1);%等距节点,节点数5y1=sin(pi.*x1);f01=Newtominter(x1,y1,x0);x2=linspace(0,1,n2);%等距节点,节点数10y2=sin(pi.*x2);f02 = Newtominter(x2,y2,x0);x3=linspace(0,1,n3);%等距节点,节点数15y3=sin(pi.*x3);f03= Newtominter(x3,y3,x0);plot(x0,y0,'-r')%原图hold onplot(x0,f01,'-g')%5个节点plot(x0,f02,'-k')%10个节点plot(x0,f03,'-b')%15个节点legend('原图','5个节点Newton插值多项式','10个节点Newton插值多项式','15个节点Newton插值多项式')运行结果：取不同的节点做牛顿插值。

相关主题