矩阵与数值分析上机实验题及程序
北航数值分析计算实习题目二 矩阵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 -==-∑，在计算每个系数时会产生相应的计算误差。
大连理工大学矩阵与数值分析上机作业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时的结果，这里就不罗列了。
矩阵与数值分析
2011级工科硕士研究生《矩阵与数值分析》课程数值实验报告教学班号：1任课老师：张宏伟姓名：xxx院系：机械工程学院学号：21104218一、 对于数列11111,,,,,392781，有如下两种生成方式1、首项为01a =，递推公式为11,1,2,3nn a a n -== ； 2、前两项为0111,3a a ==，递推公式为1210,2,3,3n n n a a a n --=-= ； 给出利用上述两种递推公式生成的序列的第50项。
1题程序: clear; clc;a=linspace(0,0,50); k=1;a(1)=1; %a0=a(1),a1=a(2)依次类推....a49=a(50)。
for k=1:49a(k+1)=1/3*a(k); endformat short e a 运行结果 a =Columns 1 through 61.0000e+000 3.3333e-001 1.1111e-001 3.7037e-002 1.2346e-002 4.1152e-003Columns 7 through 121.3717e-003 4.5725e-004 1.5242e-004 5.0805e-005 1.6935e-005 5.6450e-006Columns 13 through 181.8817e-006 6.2723e-0072.0908e-007 6.9692e-008 2.3231e-008 7.7435e-009Columns 19 through 242.5812e-009 8.6039e-010 2.8680e-010 9.5599e-0113.1866e-011 1.0622e-011Columns 25 through 303.5407e-012 1.1802e-012 3.9341e-013 1.3114e-0134.3712e-014 1.4571e-014Columns 31 through 364.8569e-015 1.6190e-0155.3966e-016 1.7989e-016 5.9962e-017 1.9987e-017Columns 37 through 426.6625e-018 2.2208e-0187.4027e-019 2.4676e-0198.2253e-020 2.7418e-020Columns 43 through 489.1392e-021 3.0464e-021 1.0155e-021 3.3849e-022 1.1283e-022 3.7610e-023Columns 49 through 501.2537e-023 4.1789e-024clear;clc;a=linspace(0,0,50);k=3;a(1)=1; %a0=a(1),a1=a(2)依次类推....a49=a(50)。
矩阵与数值分析
2013级工科硕士研究生《矩阵与数值分析》课程数值实验题目一、设622101NNjSj==-∑，分别编制从小到大和从大到小的顺序程序分别计算100001000000,S S并指出两种方法计算结果的有效位数。
Matlab程序如下：function [si,sd]=S(N)format long;si=0;sd=0;for j=N:-1:2si=1.0e6/(j^2-1)+si;endfor j=2:Nsd=1.0e6/(j^2-1)+sd;endend在matlab命令窗口中输入：[si,sd]=S(10000)运行结果：si =7.499000049995000e+005sd =7.499000049994994e+005在matlab命令窗口中输入：[si,sd]=S(1000000)运行结果：si =7.499990000005000e+005sd =7.499990000005200e+0051结果分析：si为从大到小的顺序求和的值，sd为从小到大的顺序求和的值。
当N分别为10000和1000000时，si分别为7.499000049995000e+005和7.499990000005000e+005，可以看出这两个数的有效值均为13位；而sd分别为7.499000049994994e+005和7.499990000005200e+005，这两个数的有效值均为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 --------+-⎛⎫ ⎪- ⎪= ⎪- ⎪-+⎝⎭O O O O ， 其中1211112n T --⎛⎫ ⎪- ⎪= ⎪- ⎪-⎝⎭O O O O 为1n -阶矩阵，1n I -为1n -阶单位矩阵，迭代法计算停止的条件为：101210k k x x -+-<，给出10,100,1000n =时的不同迭代步数． 求解系数矩阵A 和向量b 的Matlab 程序如下：function [A,b,x0]=jz(n)for i=1:n-1 %此处n 赋值即可，如n=100for j=1:n-1if(i==j)T(i,j)=2;endif(i==(j+1))T(i,j)=-1;endif(i==(j-1))T(i,j)=-1;endendendI=eye(n-1);k=1;for mm=1:(n-1)A(k:(k+n-2),k:(k+n-2))=T+2*I;k=k+n-1;endk=1;for xx=1:(n-2)A(k:(k+n-2),(k+n-1):(k+2*n-3))=-I;A((k+n-1):(k+2*n-3),k:(k+n-2))=-I;k=k+n-1;endb=randn((n-1)^2,1);x0=zeros((n-1)^2,1);Jacobi 迭代法Matlab 程序如下：function n=jacobi(A,b,x0)eps= 1.0e-10;M = 10000;D=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵B=D\(L+U);f=D\b;x=B*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=B*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endendGauss-Seidel迭代法Matlab程序如下：function n=gauseidel(A,b,x0)eps= 1.0e-10;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+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！');return;endend在matlab命令窗口中输入：[A,b,x0]=jz(10);J10=jacobi(A,b,x0)G10=gauseidel(A,b,x0)[A,b,x0]=jz(20);J20=jacobi(A,b,x0)G20=gauseidel(A,b,x0)[A,b,x0]=jz(30);J30=jacobi(A,b,x0)G30=gauseidel(A,b,x0)运行结果：J10 =436；G10 =214J20 =1810；G20 =913J30 =3990；G30 =2001结果分析：N=10且M=1000时，Jacobi 迭代法和Gauss —seidel 迭代法的迭代次数分别为436和214；N=20且M=1000时，Jacobi 迭代法和Gauss —seidel 迭代法的迭代次数分别为1810和913；N=30且M=1000时 ，Jacobi 和Gauss-seidel 算法的迭代次数分别为3990和2001次。
大连理工大学矩阵与数值分析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
线性方程组部分实验题11. 问题描述2. 问题分析2.1. Hilbert 矩阵Hilbert 矩阵是一种著名的“坏条件”矩阵，及病态矩阵，该矩阵任意一个元素的微小变化都能引起巨大的误差（如求逆等）。
该矩阵的元素的数学表达式是a(i ，j)=1/(i+j-1)。
这里我们用matlab 自带的产生Hilbert 矩阵的函数hilb(n)计算（n=1、 2......15）。
同时可以用MATLAB 计算。
2.2. Gauss 消去法（LU 分解）和Cholesky 分解方法高斯消去法主要有顺序消去和回代过程，是线性方程最基本的解法之一。
对n n A R ⨯∈，若0detA ≠，0(1,2,,)i i n ∆≠=⋅⋅⋅，则存在唯一的单位下三角形矩阵L 和非奇异上三角形矩阵U 使得A=LU 。
Cholesky 分解方法，也称平方根法。
对n n A R ⨯∈，假设A 对称正定，则存在唯一对角元素为正数的下三角矩阵L ，使得T A=LL 。
2.3. 条件数计算根据条件数定义：1222()n n n cond H H H -= ；12221n i i X x =⎛⎫== ⎪ ⎪⎝⎭∑ 3. 求解方案设计3.1. 线性方程组的求解及误差在这里为了比较，选择n=2、5、8、10、15 五种维数进行计算分析。
按照题意，我们先设x 为分量全1的向量，求出b ，然后将H 和b 作为已知量，求x ，与设定的参考解对比。
Gauss 消去法、Cholesky 分解法的Matlab 实现程序如下：----------------------------------------------------------------------------------------------------------%Guass 法***************************************************************************①Guass.m function guass(n)n=input('请输入系数矩阵的维数n='); %输入nH=hilb(n); %构造Hilbert 矩阵 x_exact=ones(n ，1); %构造x 列向量 b=H*x_exact; %构造b 向量 x=Doolittle(H ，b)a=input('是否继续，继续请按1，结束请按0:'); if a==1clear;clc;close all; guass(n); else end②Doolittle.mfunction x=Doolittle(A ，b) % LU 分解求Ax=b 的解 N=size(A); n=N(1);L=eye(n ，n);%L 的对角元素为1 U=zeros(n ，n); %U 为零矩阵 U(1，1:n)=A(1，1:n)%U 第一行L(1:n，1)=A(1:n，1)/U(1，1)%L第一列for k=2:nfor i=k:nU(k，i)=A(k，i)-L(k，1:(k-1))*U(1:(k-1)，i)endfor j=(k+1):nL(j，k)=(A(j，k)-L(j，1:(k-1))*U(1:(k-1)，k))/U(k，k)endendy=solveDownTriangle(L，b);%调用下三角系数矩阵求解函数x=solveUpTriangle(U，y);%调用上三角系数矩阵求解函数③solveDownTriangle.mfunction x=solveDownTriangle(A，b)%求下三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=1:nif (i>1)s=A(i，1:(i-1))*x(1:(i-1)，1);elses=0;endx(i，1)=(b(i)-s)/A(i，i);end④solveUpTriangle.mfunction x=solveUpTriangle(A，b)% 求上三角系数矩阵的线性方程组Ax=bN=size(A);n=N(1);for i=n:-1:1if (i<n)s=A(i，(i+1):n)*x((i+1):n，1);elses=0;endx(i，1)=(b(i)-s)/A(i，i);end%Cholosky分解方法******************************************************************** ---------------------------------------------------------------------------------------------------------- %利用对称正定矩阵之Cholesky分解求解线性方程组Ax=bfunction x=Chol_Solve(A，b)clear;clc% A=input('输入A=');% b=input('输入b=');n=input('请输入系数矩阵的维数n=');%定义hilbert矩阵A=hilb(n);x_exact=ones(n，1);b=A*x_exact;n=length(b);l=Cholesky(A);x=ones(n，1);y=ones(n，1);y(1)=b(1)/l(1，1);%求Ly=bfor i=2:nz=0;for k=1:i-1z=z+l(i，k)*y(k);endy(i)=(b(i)-z)/l(i，i);endx(n)=y(n)/l(n，n);%求L'x=yfor i=n-1:-1:1z=0;for k=i+1:nz=z+l(k，i)*x(k);endx(i)=(y(i)-z)/l(i，i);endxwc=x-x_exact %误差计算a=input('是否继续，继续请按1，结束请按0:');if a==1Chol_Solve(A，b);else end%对称正定矩阵之Cholesky分解法：function L=Cholesky(A)% n=input('请输入系数矩阵的维数n=');%定义hilbert矩阵% A=hilb(n);n=length(A);L=zeros(n);for k=1:ndelta=A(k，k);for j=1:k-1delta=delta-L(k，j)^2;endif delta<1e-10return;endL(k，k)=sqrt(delta);for i=k+1:nL(i，k)=A(i，k);for j=1:k-1L(i，k)=L(i，k)-L(i，j)*L(k，j);endL(i，k)=L(i，k)/L(k，k);endend3.2.条件数的计算根据2.3叙述原理，写出条件数的计算程序如下：%tiaojianshu.mm=input ('input m:=') ;N=[1:m];for i=1:mn=N(i);H=hilb(n);k=cond(H);disp('矩阵的阶数')disp(n)disp('矩阵')disp(H)disp('矩阵的条件数')disp(k)end4.结果与误差分析4.1.线性方程组计算结果分别用Gauss法和Cholesky分解法计算线性方程组的结果，并进行误差分析结果如下：4.2.条件数计算结果4.3. 结果分析取不同的n 值，得到如下结果:对于Guass 法，可以看出来，随着n 的增大，求解结果误差变大，这是因为随着n 增大，系数矩阵的条件数变大，微小的扰动就容易造成很大的误差。
2016矩阵与数值分析上机作业满分
第一题T考虑计算给定向量的数：输入向量*=（和叼…卄龙J,输出klbll刈d|x||g。
请编制一个通用程序，并用你编制的程序计算如下向量的数：,丄)_ y = (1*2n)r nI 2 3对n = 10，100，1000甚至更大的n计算其数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？代码：function Problem1()clcN=in put(' in put N=');X=zeros(1,N);Y=1:N;for i=1:NX(i)=1/i;endx1=0;x2=0;y1=0;y2=0;for i=1:Nx1=x1+abs(X(i)); x2=x2+X(i)*X(i);y仁y1+abs(Y(i)); y2=y2+Y(i)*Y(i);endx1;x2=sqrt(x2);x3=max(abs(X));y1;y2=sqrt(y2);y3=max(abs(Y));fprintf('x 的1 数是%.2f\n',x1); fprintf('x 的2 数是%.2f\n',x2);fprintf('x 的无穷数是%.2f\n',x3); fprintf('y 的1 数是%.2f\n',y1);fprintf('y 的2 数是%.2f\n',y2); fprintf('y 的无穷数是%.2f\n',y3);结果：in put N=10x的1数是2.93x的2数是1.24x的无穷数是1.00y的1数是55.00y的2数是19.62y的无穷数是10.00input N=100x的1数是5.19x的2数是1.28x的无穷数是1.00y的1数是5050.00y的2数是581.68y的无穷数是100.00in put N=1000x的1数是7.49x的2数是1.28x的无穷数是1.00y 的1 数是500500.00y 的2 数是18271.11y的无穷数是1000.00第二题==旳阿】f可考虑，其中定义f(0)=1，此时f(x)是连续函数。
- 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和拟合多项式的图形。