矩阵与数值分析上机实验题及程序
北航数值分析计算实习题目二 矩阵QR分解
数值分析实习二院（系）名称航空科学与工程学院专业名称动力工程及工程热物理学号SY0905303学生姓名解立垚1. 题目试用带双步位移QR 的分解法求矩阵A=[a ij ]10*10的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知()sin 0.50.2,1.5cos 1.2,ij i j i j a i j i j ⎧⎫+≠⎪⎪=⎨⎬+=⎪⎪⎩⎭(),1,2,...,10i j =。
说明：1、求矩阵特征值时，要求迭代的精度水平为1210ε-=。
2、打印以下内容：算法的设计方案；全部源程序(要求注明主程序和每个子程序的功能)； 矩阵A 经过拟上三角话之后所得的矩阵()1n A -；对矩阵()1n A-进行QR 分解方法结束后所得的矩阵；矩阵A 的全部特征值()(),1,2,......10i i iR I i λ=，和A 的相应于实特征值的特征向量；其中()(),.i e i m i R R I I λλ==如果i λ是实数，则令0.i I =3、采用e 型输出数据，并且至少显示12位有效数字。
2. 算法设计方案本题采用带双步位移的QR 分解方法。
为了使程序简洁，自定义类Xmatrix ，其中封装了所需要的函数方法。
在Xmatrix 类中封装了运算符重载的函数，即定义了矩阵的加、减、乘、除、数乘运算及转置运算(T())。
同时为了避免传递数组带来的额外内存开销，使用引用(&)代替值传递，以节省内存空间，避免溢出.（1）此程序的主要部分为Xmatrix 中的doubleQR()方法，具体如下：Step1：使用矩阵拟上三角化的算法将A 化为拟上三角阵A （n-1）(此处调用Xmatrix 中的preQR()方法)Step2：令121,,10k m n ε-===， 其中k 为迭代次数。
Step3：如果,1m m a ε-≤，则得到A 的一个特征值,m m a ，令1m m =-，goto Step4;否则goto Step5.Step4： 如果1m =，则得到A 的一个特征值11a ，goto Step11;如果0m =，则goto Step11;如果1m >，则goto Step3;Step5(Step6)：如果2m =，则得到A 的两个特征值12s s 和（12s s 和为右下角两阶子阵对应的特征方程21,1,()det 0m m m m a a D λλ---++=的两个根。
矩阵分析与数值分析实验报告
《矩阵分析与数值分析》实验报告院系：姓名：学号：所在班号：任课老师：一．设错误！未找到引用源。
，分别编制从小到大和从大到小的顺序程序计算错误！未找到引用源。
并指出有效位数。
程序如下：function sum3j=input('请输入求和个数 "j":');A=0;B=0;double B;double A;for n=2:jm=n^2-1;t=1./m;A=A+t;enddisp('从小到大：')s=Afor n=j:-1:2m=n^2-1;t=1./m;B=B+t;enddisp('从大到小:')s=B运行结果：>> sum3请输入求和个数 "j":100从小到大：s =0.740049504950495从大到小:s =0.740049504950495>> sum3请输入求和个数 "j":10000从小到大：s =0.749900004999506从大到小:s =0.749900004999500>> sum3请输入求和个数 "j":1000000从小到大：s =0.749999000000522从大到小:s =0.749999000000500二、解线性方程组1．分别Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组。
⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----000121001210012100124321x x x x 迭代法计算停止的条件为：6)()1(3110max -+≤≤<-k j k j j x x 。
解：（1）Jacobi 迭代法程序代码： function jacobi(A, b, N) clc;clear;A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2]; b=[-1 0 0 0]'; N=100;n = size(A,1); D = diag(diag(A)); L = tril(-A,-1); U = triu(-A,1); Tj = inv(D)*(L+U); cj = inv(D)*b; tol = 1e-06; k = 1;format longx = zeros(n,1); while k <= Nx(:,k+1) = Tj*x(:,k) + cj;disp(k); disp('x = ');disp(x(:,k+1)); if norm(x(:,k+1)-x(:,k)) < toldisp('The procedure was successful')disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1)); break endk = k+1; end结果输出The procedure was successfulCondition ||x^(k+1) - x^(k)|| < tol was met after k iterations 60 x =0.799998799067310.599998427958700.399998056850090.19999902842505（2）Gauss-Seidel迭代法程序代码：function gauss_seidel(A, b, N)clc;clear;A=[-2 1 0 0;1 -2 1 0;0 1 -2 1;0 0 1 -2];b=[-1 0 0 0]';N=100;n = size(A,1);D = diag(diag(A));L = tril(-A,-1);U = triu(-A,1);Tg = inv(D-L)*U;cg = inv(D-L)*b;tol = 1e-06;k = 1;x = zeros(n,1);while k <= Nx(:,k+1) = Tg*x(:,k) + cg;disp(k); disp('x = ');disp(x(:,k+1));if norm(x(:,k+1)-x(:,k)) < toldisp('The procedure was successful')disp('Condition ||x^(k+1) - x^(k)|| < tol was met after k iterations') disp(k); disp('x = ');disp(x(:,k+1));breakendk = k+1;end结果输出The procedure was successfulCondition ||x^(k+1) - x^(k)|| < tol was met after k iterations31x =0.799999213979350.599998971085610.399999167590770.199999583795392. 用Gauss列主元消去法、QR方法求解如下方程组：⎪⎪⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---017435222331325212214321x x x x (1)Gauss 列主元消去法 程序代码：function x=Gaussmain(A,b) clc;clear; format longA=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3]; b=[4 7 -1 0]'; N=length(A); x=zeros(N,1); y=zeros(N,1); c=0; d=0;A(:,N+1)=b; for k=1:N-1 for i=k:4if c<abs(A(i,k))d=i;c=abs(A(i,k)); end endy=A(k,:);A(k,:)=A(d,:); A(d,:)=y; for i=k+1:N c=A(i,k);for j=1:N+1A(i,j)=A(i,j)-A(k,j)*c/A(k,k); end end endb=A(:,N+1);x(N)=b(N)/A(N,N); for k=N-1:-1:1x(k)=(b(k)-A(k,k+1:N)*x(k+1:N))/A(k,k); end结果输出 ans =18.00000000000000 -9.571428571428576.00000000000000-0.42857142857143（2）QR方法：程序代码function QR(A,b)clc;clear;format longA=[1 2 1 2;2 5 3 -2;-2 -2 3 5;1 3 2 3];b=[4 7 -1 0]';[Q,R]=qr(A)y=Q'*b;x=R\y结果输出Q =-0.31622776601684 -0.04116934847963 -0.75164602800283 0.57735026918962 -0.63245553203368 -0.49403218175557 -0.15032920560056 -0.57735026918963 0.63245553203368 -0.74104827263336 -0.22549380840085 -0.00000000000000 -0.31622776601684 -0.45286283327594 0.60131682240226 0.57735026918963 R =-3.16227766016838 -6.00832755431992 -0.94868329805051 2.84604989415154 0 -2.42899156029822 -4.65213637819829 -4.15810419644272 0 0 -0.67648142520255 -0.52615221960200 0 0 0 4.04145188432738 x =17.99999999999989-9.571428571428515.99999999999997-0.42857142857143三、非线性方程的迭代解法1．用Newton迭代法求方程()06cos22x=-++=-xexf x的根，计算停止的条件为：6110-+<-kkxx；编程如下：function newton(f,df,x,a,a0)syms xf=input('please enter your equation:') a0=input('please enter you x(0):');df=diff(f)e=1e-6;a1=a0+1;N=0;while abs(a1-a0)>ea=a0-subs(f,a0)/subs(df,a0); a1=a0; a0=a; N=N+1; endfprintf('a=%0.6f',a) N运行结果： >> newtonplease enter your equation:exp(x)+2^(-x)+2*cos(x)-6 f =exp(x)+2^(-x)+2*cos(x)-6 please enter you x(0):2df =exp(x)-2^(-x)*log(2)-2*sin(x) a=1.829384 N =42．利用Newton 迭代法求多项式07951.2954.856.104.5x 234=+-+-x x x的所有实零点，注意重根的问题。
西南交通大学2018-2019数值分析Matlab上机实习题
西南交通⼤学2018-2019数值分析Matlab上机实习题数值分析2018-2019第1学期上机实习题f x，隔根第1题.给出⽜顿法求函数零点的程序。
调⽤条件：输⼊函数表达式()a b，输出结果：零点的值x和精度e，试取函数区间[,]，⽤⽜顿法计算附近的根，判断相应的收敛速度，并给出数学解释。
1.1程序代码：f=input('输⼊函数表达式：y=','s');a=input('输⼊迭代初始值：a=');delta=input('输⼊截⽌误差：delta=');f=sym(f);f_=diff(f); %求导f=inline(f);f_=inline(f_);c0=a;c=c0-f(c0)/f_(c0);n=1;while abs(c-c0)>deltac0=c;c=c0-f(c0)/f_(c0);n=n+1;enderr=abs(c-c0);yc=f(c);disp(strcat('⽤⽜顿法求得零点为',num2str(c)));disp(strcat('迭代次数为',num2str(n)));disp(strcat('精度为',num2str(err)));1.2运⾏结果：run('H:\Adocument\matlab\1⽜顿迭代法求零点\newtondiedai.m')输⼊函数表达式：y=x^4-1.4*x^3-0.48*x^2+1.408*x-0.512输⼊迭代初始值：a=1输⼊截⽌误差：delta=0.0005⽤⽜顿法求得零点为0.80072迭代次数为14精度为0.00036062⽜顿迭代法通过⼀系列的迭代操作使得到的结果不断逼近⽅程的实根，给定⼀个初值，每经过⼀次⽜顿迭代，曲线上⼀点的切线与x轴交点就会在区间[a,b]上逐步逼近于根。
上述例⼦中，通过给定初值x=1，经过14次迭代后，得到根为0.80072，精度为0.00036062。
数值分析大作业
数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。
但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。
对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。
求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。
使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。
求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。
大连理工大学矩阵与数值分析上机作业
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.设, 其精确值为.(1)编制按从⼤到⼩的顺序, 计算的通⽤程序(2)编制按从⼩到⼤的顺序, 计算的通⽤程序(3)按两种顺序分别计算并指出有效位数（编制程序时⽤单精度）(4)通过本上机题，你明⽩了什么从⼩到⼤,代码:%1---SN = %N = input('please input a number(N>=2)')if(N < 2)disp('wrong number')elseS = 0;for j = 2:1:NS = S + 1/(j^2 -1);enddisp('S:')disp(S)end结果please input a number(N>=2)10^2N =100S:7.4005e-001>> clearplease input a number(N>=2)10^4N =10000S:7.4990e-001>> clearplease input a number(N>=2)10^6N =1000000S:7.5000e-001>>从⼤到⼩代码:%1---SN = %eps('single')N = input('please input a number(N>=2)') if(N < 2) disp('wrong number')elseS = 0;for j = N:-1:2S = S + 1/(j^2 -1);enddisp('S:')disp(S)end结果please input a number(N>=2)10^2N =100S:7.4005e-001>> clearans =1.1921e-007please input a number(N>=2)10^4N =10000S:7.4990e-001>> clearans =1.1921e-007please input a number(N>=2)10^6N =1000000S:7.5000e-001(4)计算的顺序影响结果。
(完整word版)数值分析课程设计实验二
实验二2．1一、题目：用高斯消元法的消元过程作矩阵分解。
设20231812315A ⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦消元过程可将矩阵A 化为上三角矩阵U ，试求出消元过程所用的乘数21m 、31m 、31m 并以如下格式构造下三角矩阵L 和上三角矩阵U(1)(1)212223(2)313233120231,1L m U a a m m a ⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦验证：矩阵A 可以分解为L 和U 的乘积，即A =LU 。
二、算法分析：设矩阵111213212223313233a a a A a a a a a a ⎛⎫ ⎪= ⎪ ⎪⎝⎭，通过消元法可以将其化成上三角矩阵U ，具体算法如下： 第1步消元：111111(1)22112(1)331130,0;;2,3;i i i i i i i i a m a a a a m a i a a m a +=≠⎧⎪=+=⎨⎪=+⎩ 得到111213(1)(1)12223(1)(1)323300a a a A a a a a ⎛⎫ ⎪= ⎪ ⎪⎝⎭第2步消元：(1)(1)(1)32322222(2)(1)(1)333332230,0;;a m a a a a m a ⎧+=≠⎪⎨=+⎪⎩ 得到的矩阵为111213(1)(1)22223(2)33000a a a A a a a ⎛⎫ ⎪= ⎪ ⎪⎝⎭三、程序及运行结果b1.mA=[20 2 3;1 8 1;2 -3 15];for i=1:2M(i)=A(i+1,1)/A(1,1);endfor j=2:3A1(j,2)=A(j,2)-M(j-1)*A(1,2);A1(j,3)=A(j,3)-M(j-1)*A(1,3);endM(3)=A1(3,2)/A1(2,2);A1(3,2)=0;A1(3,3)=A1(3,3)-M(3)*A1(2,3);M,A1运行结果为：M =0.0500 0.1000 -0.4051A1 =0 0 00 7.9000 0.85000 0 15.0443所以：10020230.051007.90.850.10.405110015.0443L U ⎛⎫⎛⎫ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪-⎝⎭⎝⎭验证：L=[1 0 0;0.05 1 0;0.1 -0.4051 1];U=[20 2 3;0 7.9 0.85;0 0 15.0443];A1=L*UA1 =20.0000 2.0000 3.00001.0000 8.0000 1.00002.0000 -3.0003 15.0000四、精度分析因为根据LU 的递推公式可知，L ，U 分别为下三角和上三角矩阵，其中L 不在对角线上的元素值为111()k ik ik is sk s kk l a l u u -==-∑，在计算每个系数时会产生相应的计算误差。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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消去法计算结果严重偏离准确值成为错解，列主元消去法依然能得到正确解。
这是因为Gauss消去法中有小主元做除数，在计算过程中的舍入误差会对解产生极大影响，从而导致错误。
列主元消去法则避免了这种情况。
2. 设有方程组Ax b =，其中2020A R ⨯∈，30.50.250.530.50.250.50.50.250.530.50.250.53A --⎛⎫⎪-- ⎪ ⎪-- ⎪= ⎪ ⎪-- ⎪-- ⎪ ⎪--⎝⎭(a) 选取不同的初始向量(0)x和不同的右端向量b ，给定迭代误差要求，用Jacobi 和Gauss-Seidel 迭代法计算，观测得出的迭代向量序列是否收敛。
若收敛，记录迭代次数，分析计算结果并得出你的结论。
(b) 选定初始向量初始向量(0)x和不同的右端向量b ，如取(0)0,,(1,1,1)T xb Ae e === 。
将A 的主对角线元素成倍增长若干次，非主对角元素不变，每次用Jacobi 法计算，要求迭代误差满足(1)()610k k xx +-∞-<，比较收敛速度，分析现象并得出你的结论。
(a)Matlab 编程： Jacobi 迭代法：function x=Jacobi(b,x0) A=zeros(20); for i=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25; endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5; D=diag(diag(A));L=tril(A-D);U=triu(A-D); B=D\(L+U);f=D\b; x=x0;i=0;ep=10^-5; while i<1000 y=x;x=B*x+f;i=i+1; if norm(x-y,inf)<ep break; end endif i<1000disp('迭代收敛且迭代次数为') i elsedisp('迭代不收敛')endendGauss-Seidel迭代法：function x=G_S(b,x0)A=zeros(20);for i=1:18A(i,i)=3;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3;A(20,20)=3;A(19,20)=-0.5;A(20,19)=-0.5;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=(D-L)\U;f=(D-L)\b;x=x0;i=0;ep=10^-5;while i<1000y=x;x=B*x+f;i=i+1;if norm(x-y,inf)<epbreak;endendif i<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endend计算结果：x和不同的右端向量b：设误差要求为1×10^-5任选了两组不同的初始向量(0)b=[1,1,1,1,1,1,1,1,2,3,4,1,4,2,3,1,3,1,4,6]' x0=[2,3,1,4,4,5,5,3,1,3,1,1,1,1,1,1,1,1,1,1]' b=[1,4,5,12,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]' x0=[0,1,2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]' 在matlab内运行两个程序得到：第一组：>> Jacobi(b,x0)迭代收敛且迭代次数为i =19>> G_S(b,x0)迭代收敛且迭代次数为i =9第二组：>> Jacobi(b,x0)迭代收敛且迭代次数为i =18>> G_S(b,x0)迭代收敛且迭代次数为i =9分析：下Gauss-Seidel迭代法要比Jacobi迭代法收敛速度快。
因为Gauss-Seidel 迭代法对已经前面计算出来的信息进行了充分利用。
(b)Matlab编程：function x=Jacobi2(n)A=zeros(20);x=diag(A);e=diag(eye(20));for i=1:18A(i,i)=3*n;A(i+1,i)=-0.5;A(i,i+1)=-0.5;A(i+2,i)=-0.25;A(i,i+2)=-0.25;endA(19,19)=3*n;A(20,20)=3*n;A(19,20)=-0.5;A(20,19)=-0.5;b=A*e;D=diag(diag(A));L=tril(A-D);U=triu(A-D);B=D\(L+U);f=D\b;i=0;ep=10^-6;while i<1000y=x;x=B*x+f;i=i+1;if norm(x-y,inf)<epbreak;endendif i<1000disp('迭代收敛且迭代次数为')ielsedisp('迭代不收敛')endend计算结果：n代表主元增长倍数，选取不同的n得：>> Jacobi2(1)迭代收敛且迭代次数为i =20>> Jacobi2(2)迭代收敛且迭代次数为i =11>> Jacobi2(3)迭代收敛且迭代次数为i =9>> Jacobi2(4)迭代收敛且迭代次数为i =8>> Jacobi2(10)迭代收敛且迭代次数为i =6分析：矩阵主对角元增长倍数的变大，其收敛速度变快。
3.用迭代法求方程32⨯。
0.510-310+-=的全部根，要求误差限为8x x首先经过简单计算并结合图形得知该方程的三个单根区间为[-3,-1],[-1,0],[0,1],然后利用二分法求解。
Matlab编程：function x=f(a,b)fa=a^3+3*a^2-1;fb=b^3+3*b^2-1;if fa*fb>0disp('区间不是单根区间')return;endep=0.5*10^-8;while abs(a-b)>epc=(a+b)/2;f1=a^3+3*a^2-1; f2=c^3+3*c^2-1; j=f1*f2; if j==0 x=c; break; elseif j<0 b=c; else a=c; end x=a; end enddisp('该区间内根为') end计算结果： >> f(-3,-1) 该区间内根为 ans =-2.879385244101286>> f(-1,0)该区间内根为 ans =-0.652703646570444>> f(0,1)该区间内根为 ans =0.5320888832211494. 给定数据如下表：编制程序求三次样条插值函数在插值中点的样条函数插值，并作点集{,}i i x y 和样条插值函数的图形，满足的边界条件为（1）(0)0.8,(10)0.2.s s ''== （2）(0)(10)0.s s ''''==(1)Matlab编程：clear allx=[0 1 2 3 4 5 6 7 8 9 10];f1=0.8;f2=0.2;y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n-2);yn=[];xn=[];for k=2:nhk=x(n)-x(n-1);h=[h hk];endfor k=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[g gk];t=[t,tk];u=[u uk];endfor k=2:n-2A(k-1,k)=u(k-1);A(k,k-1)=t(k);endg(1)=g(1)-t(1)*f1;g(n-2)=g(n-2)-u(n-2)*f2;m=A\g';m=m';m=[f1 m f2];for k=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2; a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2; b2=y(k+1)/h(k)^2;c1=[1 -x(k)];c2=[1 -x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0 sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfor k=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[yn yi];xn=[xn xi];endxn=[xn x(n)];yn=[yn polyval(s(n-1,:),x(n))];plot(xn,yn, 'x');title('3次样条函数插值结果')计算结果：在Matlab命令框内输入上述程序得到插值中点的样条函数值Y和样条插值函数的图形：Y =0.398564254606255Y =1.168428726968726Y =1.871470837518841Y =2.478187922956004Y =2.873277470657021Y =3.213702194414573Y =3.084413751685133Y =2.919892798844714Y =3.149765052935493Y =3.222296989411632(2)Matlab编程：clear allx=[0 1 2 3 4 5 6 7 8 9 10];f1=0;f2=0;y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29];n=length(x);h=[];t=[];g=[];u=[];s=[];A=2*eye(n);yn=[];xn=[];for k=2:nhk=x(n)-x(n-1);h=[h hk];endfor k=2:length(h)tk=h(k)/(h(k)+h(k-1));uk=h(k-1)/(h(k)+h(k-1));gk=3*(uk*(y(k+1)-y(k))/h(k)+tk*(y(k)-y(k-1))/h(k-1));g=[g gk];t=[t,tk];u=[u uk];endfor k=1:n-2A(k+1,k+2)=u(k);A(k+1,k)=t(k);endA(1,2)=1;A(n,n-1)=1;g0=3*(y(2)-y(1))/h(1)-h(1)*f1/2;gn=3*(y(n)-y(n-1))/h(n-2)+h(n-2)*f2/2;g=[g0,g,gn];m=A\g';m=m';for k=1:n-1a1=2*y(k)/h(k)^3;a2=-2*y(k+1)/h(k)^3;a3=m(k)/h(k)^2; a4=m(k+1)/h(k)^2;b1=y(k)/h(k)^2; b2=y(k+1)/h(k)^2;c1=[1 -x(k)];c2=[1 -x(k+1)];c3=conv(c1,c1);c4=conv(c2,c2);sk1=b1*c4+b2*c3;sk2=(a1+a3)*conv(c1,c4)+(a2+a4)*conv(c2,c3);sk=[0 sk1]+sk2;s=[s;sk];Pn=poly2str(sk,'x');Y=polyval(s(k,:),(x(k)+x(k+1))/2)endfor k=1:n-1xi=[x(k):0.02:x(k+1)-0.02];yi=polyval(s(k,:),xi);yn=[yn yi];xn=[xn xi];endxn=[xn x(n)];yn=[yn polyval(s(n-1,:),x(n))];plot(xn,yn, 'x');title('3次样条函数插值结果')计算结果：在命令框内输入上述程序得到：Y =0.398428148708663Y =1.168465553874013Y =1.871459635795274Y =2.478195902944736Y =2.873256752425371Y =3.213777087353492Y =3.084134898160016Y =2.920933320005332Y =3.145881821815571Y =3.236789392730741x y的图形再次输入plot(x,y, 'x');可以得到点集{,}i i5.对下列数据作三次多项式拟合，取权系数1w ，给出拟合多项式的系数、平方误差，ix y和拟合多项式的图形。