大连理工大学矩阵与数值分析上机作业
大连理工大学矩阵分析网上作业
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; }
大连理工大学矩阵与数值分析上机作业代码
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
《矩阵与数值分析》上机大作业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综上，高斯列主元消去法可以避免小数作除数带来的误差，获得满意的数值解。
大连理工大学矩阵与数值分析上机作业
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 为输入参数
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 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插值多项式')运行结果：取不同的节点做牛顿插值。