矩阵与数值分析实习题2018秋

合集下载

数值分析课程实验设计——数值线性代数实习题

数值分析课程实验设计——数值线性代数
实习题
1. 实验目的
本实验的主要目的是进一步加深对数值线性代数的理解，熟悉
常见矩阵分解方法，并在此基础上解决实际问题。

2. 实验内容
本次实验将任务分为两个部分，分别是矩阵分解与求解线性方
程组。

2.1 矩阵分解
首先，我们需要熟悉三种常见的矩阵分解：QR分解、LU分解
和奇异值分解。

我们需要通过Python语言实现这三种分解方法，
并利用这些方法解决实际问题。

2.2 求解线性方程组
其次，我们需要学会用矩阵分解的方法来求解线性方程组。

我
们将通过两个例子来进行说明，并利用Python语言实现这些方法。

3. 实验要求
本次实验要求熟悉矩阵分解的基本方法，在此基础上解决实际问题；能够运用多种方法来求解线性方程组，并分析比较它们的优缺点。

4. 实验总结
本次实验通过矩阵分解和求解线性方程组两个部分的学习，巩固了我们对于数值线性代数的知识，并在实际问题的解决中得到了应用。

感谢老师的指导，我们会在今后的学习中持续探索数值分析方面的知识。

北航数值分析计算实习题目二矩阵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.给定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消去法计算结果严重偏离准确值成为错解，列主元消去法依然能得到正确解。

矩阵与数值分析试题

1）1）已知近似值246.00a =有5位有效数字，则a 的绝对误差界为，a 的相对误差界为；2）于0,2π⎡⎤⎢⎥⎣⎦，用y=a+bx 做()sin f x x =最佳平方逼近，则法方程组为：；3）设71057⎛⎫= ⎪⎝⎭A ，1=A ，()1cond =A ；4）为了减少运算次数，应将表达式.425432168116171814131x x x x x x x x ++--+---改写为_ ______；5）已知(0)1,(1)3,(2)5,f f f ===则均差[0,1,2]f = ，对应于x 0=0插值基函数()0l x = ；6）此数值求积公式2110116x e dx e --⎛⎫≈ ⎪⎝⎭⎰的代数精度为：；7) 求解1u u t e -'=-+-的隐式Euler 公式：； 8) 用二分法求方程3()2510f x x x =--=在区间[1,3]内的根，进行一步后根所在区间为___ ___。

9)1225⎛⎫=⎪⎝⎭A 的TLL 分解为：； 10) [0,1]上以1()lnx xρ=权函数的正交多项式()0x φ= ，()1x φ= 。

11）0=x 是()10x f x x e =--=的根，则具有平方收敛的迭代公式为：。

12）将向量=x 221⎛⎫⎪⎪ ⎪⎝⎭变换为向量=y ⎪⎪⎪⎭⎫ ⎝⎛003的正交矩阵H 为；姓名：学号：院系：级班装订线二、计算题1．（15分）如下求解初值问题00)(),,(u t u u t f u =='的线性二步法21(3)2n n n n hu u f f ++=++①确定出它的阶p 、局部截断误差主项和收敛性，求出其绝对稳定区间； ②给出上述方法求解方程：40u u '=-，1)0(=u ，的步长h 的取值范围。

2．（15分）确定0x ，0A ，1x ，1A 使得求积公式()()()1200111x f x dx A f x A f x -≈+⎰的代数精度m 达到最高，试问m 是多少？取()2xf x e -=，利用所求得的公式计算出数值解。

数值分析计算实习题

数值分析计算实习题-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN《数值分析》计算实习题姓名：学号：班级：第二章1、程序代码Clear;clc;x1=[ ];y1=[ ];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i-1)*df(i);endP4=vpa(sum(d),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=;q1=q(1,:)*[^3;^2;;1];q1=vpa(collect(q1),5)q2=q(1,:)*[^3;^2;;1];q2=vpa(collect(q2),5)q3=q(1,:)*[^3;^2;;1];q3=vpa(collect(q3),5)q4=q(1,:)*[^3;^2;;1];q4=vpa(collect(q4),5)%求解并化简多项式2、运行结果P4 =*x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - + q1 =- *x^3 + *x^2 - *x +q2 =- *x^3 + *x^2 - *x + q3 =- *x^3 + *x^2 - *x + q4 =- *x^3 + *x^2 - *x +3、问题结果4次牛顿差值多项式4()P x = *x - *(x - *(x - - *(x - *(x - *(x - - *(x - *(x - *(x - *(x - +三次样条差值多项式()Q x0.10.20.30.40.50.60.70.80.910.40.50.60.70.80.911.1323232321.33930.803570.40714 1.04,[0.2,0.4]1.3393 1.60710.88929 1.1643,[0.4,0.6]1.3393 2.4107 1.6929 1.4171,[0.6,0.8]1.3393 3.21432.8179 1.8629,[0.8,1.0]x x x x x x x x x x x x x x x x ⎧-+-+∈⎪-+-+∈⎪⎨-+-+∈⎪⎪-+-+∈⎩第三章1、程序代码Clear;clc; x=[0 1]; y=[1 ];p1=polyfit(x,y,3)%三次多项式拟合 p2=polyfit(x,y,4)%四次多项式拟合 y1=polyval(p1,x);y2=polyval(p2,x);%多项式求值plot(x,y,'c--',x,y1,'r:',x,y2,'y-.')p3=polyfit(x,y,2)%观察图像，类似抛物线，故用二次多项式拟合。

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

传播优3^ Word版文档•希望对您有帮助.可双击去除!矩阵与数值分析上机作业学校：大连理工大学学院：班级：姓名：学号：授课老师：注：编程语言Mat lab1.考虑计算给定向量的范数：输入向量x = (ri,x2,•■- ,-r n)r»输出胡I” ||创2，||広||oc・请编制一个通用程序，并用你编制的程序计算如下向量的范数：对“ =10, 100, 1000甚至更大的八计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？程序：Norm, m函数function s二Norm(x,m)给求向量x的范数師取1,2, in f分别表示仁2,无穷范数n=length(x);s=0;switch mcase 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的范数Testi.mclear all;clc;n1=10；n2=100；n3=1000;x1=1./[1:n1],;x2=1./[1：n2]t ;x3=1./[1:n3]f;y1 = [1:n1]t;y2=[1:n2]-;y3=[1:n3],;disp(*n=10 时J ;dispC x 的 1 一范数：*) ;disp(Norm(x1, 1)); dispC x 的 2-范数：*) ;disp(Norm(x1,2)); dispC x 的无穷-范数:*) ;disp(Norm(x1, inf)); dispCy 的1-范数：’);disp(Norm(y1, 1)); dispC'y 的 2-范数:*) ;disp(Norm(y1,2));disp(* y 的无穷-范数：*) ;disp(Norm(y1, inf)); disp(*n=100 时J;dispC x 的1-范数：*) ;disp(Norm(x2, 1));dispCx 的 2-范数:*) ;disp(Norm(x2, 2)); disp(' x 的无穷-范數:*) ;disp(Norm(x2, inf)); disp(' y 的 1 一范数：*) ;disp(Norm(y2, 1)); dispC y 的 2-范数：*) ;disp(Norm(y2, 2));disp(* y 的无穷-范数:*) ;disp(Norm(y2, inf)); dispCn=1000 时’);disp(' x 的 1-范数：*) ;disp(Norm(x3, 1)); dispC x 的 2-范数：*) ;disp(Norm(x3, 2));dispC x 的无穷-范數:*) ;disp(Norm(x3, inf)); dispC y 的 1-范数：*) ;disp(Norm(y3, 1));dispC y 的 2-范数:*) ;disp(Norm(y3, 2)); dispC y 的无穷-范数:*) ;disp(Norm(y3, inf));y 的1-范数:500500： y 的2-范数：1・8272+004； y 的无穷-范数:10002. 考虑y = /(x)=芈也，其中定义/(O) = b 此时/仗)是连续函数.用此＞公式计算当10-込10-】5］时的函数值，画出图像。

数值分析计算实习题答案

数值分析计算实习题答案数值分析计算实习题答案数值分析是一门研究如何利用计算机对数学问题进行近似求解的学科。

在数值分析的学习过程中，实习题是一种重要的学习方式，通过实践来巩固理论知识，并培养解决实际问题的能力。

本文将为大家提供一些数值分析计算实习题的答案，希望能够帮助大家更好地理解和掌握数值分析的相关知识。

一、插值与拟合1. 已知一组数据点，要求通过这些数据点构造一个一次插值多项式，并求出在某一特定点的函数值。

答案：首先，我们可以根据给定的数据点构造一个一次插值多项式。

假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个一次多项式p(x) = a0 + a1x，其中a0和a1为待定系数。

根据插值条件，我们有p(x0) = y0，p(x1) = y1。

将这两个条件代入多项式中，可以得到一个方程组，通过求解这个方程组，我们就可以确定a0和a1的值。

最后，将求得的多项式代入到某一特定点，就可以得到该点的函数值。

2. 已知一组数据点，要求通过这些数据点进行最小二乘拟合，并求出拟合曲线的表达式。

答案：最小二乘拟合是一种通过最小化误差平方和来找到最佳拟合曲线的方法。

假设给定的数据点为(x0, y0), (x1, y1)，我们可以构造一个拟合曲线的表达式y =a0 + a1x + a2x^2 + ... + anx^n，其中a0, a1, ..., an为待定系数。

根据最小二乘拟合原理，我们需要最小化误差平方和E = Σ(yi - f(xi))^2，其中yi为实际数据点的y值，f(xi)为拟合曲线在xi处的函数值。

通过求解这个最小化问题，我们就可以确定拟合曲线的表达式。

二、数值积分1. 已知一个函数的表达式，要求通过数值积分的方法计算函数在某一区间上的定积分值。

答案：数值积分是一种通过将定积分转化为数值求和来近似计算的方法。

假设给定的函数表达式为f(x)，我们可以将定积分∫[a, b]f(x)dx近似为Σwi * f(xi)，其中wi为权重系数，xi为待定节点。

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

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

北航数值分析-实习作业1(C语言详细注释)

《数值分析》计算实习作业《一》北航第一题设有501501⨯的矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=501500499321a bc b a b cc b a b ccb a bc c b a b c b a A其中.064.0,16.0);501,2,1(64.0)2.0sin()024.064.1(1.0-===--=c b i e i i a i i 矩阵的特征值)501,,2,1( =i i λ满足||min ||,501150121i i s λλλλλ≤≤=<<<试求1. 5011，λλ和s λ的值2. 的与数4015011λλκλμ-+=k 最接近的特征值)39,,2,1( =K κλi3. 的（谱范数）条件数2)A (cond 和行列式A det 要求1. 算法的设计方案（A 的所有零元素都不能存储）2. 全部源程序（详细注释）。

变量为double ，精度-1210=ε，输出为e 型12位有效数字3. 特征值s 5011，，λλλ和)39,,2,1( =K κλi 以及A cond det ,)A (2的值4. 讨论迭代初始向量的选取对计算结果的影响，并说明原因解答：1. 算法设计对于s λ满足||min ||5011i i s λλ≤≤=，所以s λ是按模最小的特征值，直接运用反幂法可求得。

对于5011，λλ，一个是最大的特征值，一个是最小的特征值，不能确定两者的绝对值是否相等，因此必须首先假设||||5011λλ≠，然后运用幂法，看能否求得一个特征值，如果可以求得一个，证明A 是收敛的，求得的结果是正确的，然后对A 进行带原点平移的幂法，偏移量是前面求得的特征值，可以求得另一个特征值，最后比较这两个特征值，较大的特征值是501λ，较小的特征值就是1λ。

如果在假设的前提下，无法运用幂法求得按模最大的特征值，即此时A 不收敛，则需要将A 进行带原点平移的幂法，平移量可以选取1，再重复上述步骤即可求得两个特征值。

矩阵与数值分析上机实习题汇总

矩阵与数值分析上机实习题汇总矩阵与数值分析上机实习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)计算的顺序影响结果。

大连理工大学矩阵与数值分析大作业题目

2014级工科硕士研究生《矩阵与数值分析》课程数值实验题目1. 方程在x=3.0附近有根，试写出其三种不同的等价形式以构成两种不同的迭代格式，再用这两种迭代求根，并绘制误差下降曲线，观察这两种迭代是否收敛及收敛的快慢2. 用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分，要求绝对误差为，并将计算结果与精确解进行比较：（1）（2）3. 使用带选主元的分解法求解线性方程组,其中，，当时．对于的情况分别求解．精确解为．对得到的结果与精确解的差异进行解释．4. 用4阶Runge-kutta 法求解微分方程t t t te e t u u u u u 222101)(,101)0(,2---+==-=' (1) 令1.0=h ，使用上述程序执行20步，然后令05.0=h ，使用上述程序执行40步(2) 比较两个近似解与精确解(3) 当h 减半时，(1)中的最终全局误差是否和预期相符？(4) 在同一坐标系上画出两个近似解与精确解．（提示输出矩阵R 包含近似解的x 和y 坐标，用命令plot(R(:,1),R(:,2))画出相应图形．）5. 设为阶的三对角方阵，是一个阶的对称正定矩阵其中为阶单位矩阵。

设为线性方程组的真解，右边的向量由这个真解给出。

(1) 用Cholesky 分解法求解该方程.(2) 用Jacobi 迭代法和Gauss-Seidel 迭代法求解该方程组，误差设为. 其中取值为4，5，6.6. 设考虑空间的一个等距划分，分点为设为插值于这些等分点上的Lagrange插值多项式。

(1)选择不断增大的分点数目画出原函数与插值多项式在的图像，并比较分析实验结果。

(2)选择重复上述的实验看其结果如何实验须知：（1）所有的数值实验的题目要求用C语言或Matlab编程；（2）实验报告内容应包括问题、程序、计算结果及分析等；（3）考试前提交实验报告；（4）本次实验成绩将占总成绩的10%。

（5）报告上要注明：所在教学班号、任课老师的姓名；报告人所在院系、学号。

数值分析计算实习第一题

直接用定义: ��(��)2 = ‖��‖2‖��−1‖2
求 A 的条件数很繁琐,需要先进行化简：
首先：
由于 A 是对称矩阵，
‖��‖2 = ��max(��)
说明：
1. 在所用的算法中，凡是要给出精度水平的ε，都取 ��=10−12。
2. 选择算法的时候应使矩阵 A 的所有零元素都不存储。
3. 打印以下内容：
(1)算法设计方案和思路。
(2)全部源程序。
(3)特征值��1,��501,��,��(��=1,2,⋯,39)以及��(��)2, det��的值(采用 e 型输出实型数，并至少显示 12 位有效数字)。
λi[16] -2.533970311130E+00 λi[38] 8.648666065193E+00
λi[17] -2.003230769563E+00 λi[39] 9.254200344575E+00
λi[18] -1.503557611227E+00 cond(A)2 1.925204273903E+03
λi[19] -9.935586060080E-01 det(A) 2.772786141752E+118

北航数值分析计算实习题目一幂法反幂法求矩阵特征值

《数值分析》计算实习题目第一题：1. 算法设计方案（1）1λ，501λ和s λ的值。

1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。

2)使用反幂法求λs ，其中需要解线性方程组。

因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。

（2）与140k λλμλ-5011=+k 最接近的特征值λik 。

通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。

（3）2cond(A)和det A 。

1）1=nλλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。

2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。

由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。

2.全部源程序#include <stdio.h>#include <math.h>void init_a();//初始化Adouble get_an_element(int,int);//取A 中的元素函数double powermethod(double);//原点平移的幂法double inversepowermethod(double);//原点平移的反幂法int presolve(double);//三角LU 分解int solve(double [],double []);//解方程组int max(int,int);int min(int,int);double (*u)[502]=new double[502][502];//上三角U 数组double (*l)[502]=new double[502][502];//单位下三角L 数组double a[6][502];//矩阵Aint main(){int i,k;double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;double lambda[40];init_a();//初始化Alambdat1=powermethod(0);lambdat2=powermethod(lambdat1);lambda1=lambdat1<lambdat2?lambdat1:lambdat2;lambda501=lambdat1>lambdat2?lambdat1:lambdat2;presolve(0);lambdas=inversepowermethod(0);det=1;for(i=1;i<=501;i++)det=det*u[i][i];for (k=1;k<=39;k++){mu[k]=lambda1+k*(lambda501-lambda1)/40;presolve(mu[k]);lambda[k]=inversepowermethod(mu[k]);}printf("------------所有特征值如下------------\n");printf("λ=%1.11e λ=%1.11e\n",lambda1,lambda501);printf("λs=%1.11e\n",lambdas);printf("cond(A)=%1.11e\n",fabs(lambdat1/lambdas));printf("detA=%1.11e \n",det);for (k=1;k<=39;k++){printf("λi%d=%1.11e ",k,lambda[k]);if(k % 3==0) printf("\n");} delete []u;delete []l;//释放堆内存return 0;}void init_a()//初始化A{int i;for (i=3;i<=501;i++) a[1][i]=a[5][502-i]=-0.064;for (i=2;i<=501;i++) a[2][i]=a[4][502-i]=0.16;for (i=1;i<=501;i++) a[3][i]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i); }double get_an_element(int i,int j)//从A中节省存储量的提取元素方法{if (fabs(i-j)<=2) return a[i-j+3][j];else return 0;}double powermethod(double offset)//幂法{int i,x1;double u[502],y[502];double beta=0,prebeta=-1000,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0;//设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;for (x1=1;x1<=501;x1++){u[x1]=0;for (int x2=1;x2<=501;x2++)u[x1]=u[x1]+((x1==x2)?(get_an_element(x1,x2)-offset):get_an_element(x1,x2))*y[x2] ;}prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);}double inversepowermethod(double offset)//反幂法{int i;double u[502],y[502];double beta=0,prebeta=0,yita=0;for (i=1;i<=501;i++)u[i]=1,y[i]=0; //设置初始向量u[]for (int k=1;k<=10000;k++){yita=0;for (i=1;i<=501;i++) yita=sqrt(yita*yita+u[i]*u[i]);for (i=1;i<=501;i++) y[i]=u[i]/yita;solve(u,y);prebeta=beta;beta=0;for (i=1;i<=501;i++) beta=beta+ y[i]*u[i];beta=1/beta;if (fabs((prebeta-beta)/beta)<=1e-12) {printf("offset=%f lambda=%f err=%e k=%d\n",offset,(beta+offset),fabs((prebeta-beta)/beta),k);break;};//输出中间过程,包括偏移量,误差,迭代次数}return (beta+offset);int presolve(double offset)//三角LU分解{int i,k,j,t;double sum;for (k=1;k<=501;k++)for (j=1;j<=501;j++){u[k][j]=l[k][j]=0;if (k==j) l[k][j]=1;} //初始化LU矩阵for (k=1;k<=501;k++){for (j=k;j<=min(k+2,501);j++){sum=0;for (t=max(1,max(k-2,j-2)) ; t<=(k-1) ; t++)sum=sum+l[k][t]*u[t][j];u[k][j]=((k==j)?(get_an_element(k,j)-offset):get_an_element(k,j))-sum;}if (k==501) continue;for (i=k+1;i<=min(k+2,501);i++){sum=0;for (t=max(1,max(i-2,k-2));t<=(k-1);t++)sum=sum+l[i][t]*u[t][k];l[i][k]=(((i==k)?(get_an_element(i,k)-offset):get_an_element(i,k))-sum)/u[k][k];}}return 0;}int solve(double x[],double b[])//解方程组{int i,t;double y[502];double sum;y[1]=b[1];for (i=2;i<=501;i++){sum=0;for (t=max(1,i-2);t<=i-1;t++)sum=sum+l[i][t]*y[t];y[i]=b[i]-sum;}x[501]=y[501]/u[501][501];for (i=500;i>=1;i--){sum=0;for (t=i+1;t<=min(i+2,501);t++)sum=sum+u[i][t]*x[t];x[i]=(y[i]-sum)/u[i][i];}return 0;}int max(int x,int y){return (x>y?x:y);}int min(int x,int y){return (x<y?x:y);}3.计算结果结果如下图所示:部分中间结果：给出了偏移量(offset)，误差(err)，迭代次数(k)4.讨论迭代初始向量的选取对计算结果的影响,并说明原因使用u[i]=1（i=1,2,...,501）作为初始向量进行迭代，可得出以上结果。

北航数值分析实习题目第一题

《数值分析B》大作业一ZY1515105 樊雪松一．算法设计方案：1.矩阵A的存储与检索将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] 。

在数组MatrixC[5][501]中检索A的带内元素a ij的方法是：A的带内元素a ij=C中的元素c i-j+2,j。

2.求解λ1，λ501，λs1、首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。

λmin即为λs；如果λ max>0,则λ501=λmax；如果λmax<0,则λ1=λmax。

2、使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ’max，如果λ max>0,则λ1=λ’max+p；如果λmax<0,则λ501=λ’max+p。

3、求解A的与数μk=λ1+k（λ501-λ1）/40 的最接近的特征值λik （k=1,2，…，39）。

使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λ ik。

4、求解A的（谱范数）条件数cond（A）2和行列式detA。

cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。

求解矩阵A的行列式，可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。

二．源程序#include<stdio.h>#include<math.h>#include<conio.h>//定义A中元素double C[5][501];double a[501];double b;double c;//声明所有函数void YaSuoJZ(double C[5][501],double a[501],double b,double c) ;//压缩矩阵函数double mifa(double C[5][501]); //幂法函数void daizhuangLU(double A[5][501]); //带状矩阵的LU分解double fanmifa(double C[5][501]);//反幂法函数//最值函数int max2(int x,int y);int max3(int x,int y,int z);int min(int x,int y);//最值函数int max2(int x,int y) //求2个数的最大值{int z;z=x>y?x:y;return(z);}int max3(int x,int y,int z) //求3个数的最大值{int w;w = z > max2(x,y)? z:max2(x,y);return(w);}int min(int x,int y) //求2个数的最小值{int z;z=x>y?y:x;return(z);}//将矩阵A压缩存储在矩阵C中void YaSuoJZ(double C[5][501],double a[501],double b,double c) {int i;for(i=0;i<=500;i++){if(i>=2) C[0][i]=c;else C[0][i]=0;if(i>=1) C[1][i]=b;else C[1][i]=0;if(i<=499) C[3][i]=b;else C[3][i]=0;if(i<=498) C[4][i]=c;else C[4][i]=0;C[2][i]=a[i];}}//幂法函数：用幂法求矩阵模最大的特征值double mifa(double C[5][501]){double u[501];double y[501]={0},η=0;double β,βk=0;double ε=1;// ε为精度double sumu=0,sumAY=0;int i,j,k=1; //k为循环次数for (i=0;i<=500;i++) //取任一非零向量u0u[i] = 1.0;while(ε>=1e-12){for(i=0;i<=500;i++) //求u(k-1)的2范数ηsumu=sumu+u[i]*u[i];η=sqrt(sumu);sumu=0;for(i=0;i<=500;i++) //求y(k-1)y[i]=u[i]/η;for(i=0;i<=500;i++) //求u(k)的各分量u[i]{for(j=max2(0,i-2);j<=min(i+2,500);j++)sumAY=sumAY+C[i-j+2][j]*y[j];u[i]=sumAY;sumAY=0;}//求幂法中的βkβ=βk; //将β(k-1)放在β中βk=0;for(i=0;i<=500;i++) //求βkβk=βk+y[i]*u[i];if(k>=2)ε=fabs(βk-β)/fabs(βk);k++;}return(βk);}//带状矩阵的LU分解void daizhuangLU(double A[5][501]){int i,j,k,m,t;double sumukj=0,sumlik=0;for(k=0;k<=500;k++){for(j=k;j<=min(k+2,500);j++) //求ukj并存在A[k-j+2][j]中{for(t=max3(0,k-2,j-2);t<=k-1;t++)sumukj=sumukj+A[k-t+2][t]*A[t-j+2][j];A[k-j+2][j]=A[k-j+2][j]-sumukj;sumukj=0;}if(k<500)for(i=k+1;i<=min(k+2,500);i++) //求lik并存在A[i-k+2][k]中{for(m=max3(0,i-2,k-2);m<=k-1;m++)sumlik=sumlik+A[i-m+2][m]*A[m-k+2][k];A[i-k+2][k]=(A[i-k+2][k]-sumlik)/A[2][k];sumlik=0;}}}//反幂法函数：用反幂法求矩阵的模最小的特征值double fanmifa(double M[5][501]){double u[501];double y[501]={0},x[501],η=0;double fβ,fβk=0;double ε=1;double fsumu=0,sumLX=0,sumUu=0;int i,t,m,k=1;for(i=0;i<=500;i++) //任取一非零向量u0u[i]=1;daizhuangLU(M); //对A进行LU分解A=LU，Au(k)=y(k-1)等价于Uu(k)=x和Lx=y(k-1) while(ε>=1e-12){for(i=0;i<=500;i++) //求u(k-1)的2范数ηfsumu=fsumu+u[i]*u[i];η=sqrt(fsumu);fsumu=0;for(i=0;i<=500;i++) //求y(k-1)y[i]=u[i]/η;for(i=0;i<=500;i++) //求中间向量xx[i]=y[i];for(i=1;i<=500;i++){for(t=max2(0,i-2);t<=i-1;t++)sumLX=sumLX+M[i-t+2][t]*x[t];x[i]=x[i]-sumLX;sumLX=0;}u[500]=x[500]/C[2][500]; //求u(k)的各分量u[i]for(i=499;i>=0;i--){for(m=i+1;m<=min(i+2,500);m++)sumUu=sumUu+M[i-m+2][m]*u[m];u[i]=(x[i]-sumUu)/M[2][i];sumUu=0;}//求反幂法中的βkfβ=fβk; //将fβ(k-1)放在fβ中fβk=0;for(i=0;i<=500;i++) //求fβkfβk=fβk+y[i]*u[i];if(k>=2)ε=fabs(1/fβk-1/fβ)/fabs(1/fβk);k++;}return(1/fβk);}//主函数void main(){int i,j,k;double λ1,λ501,λm,λm1,λm2,λs,λ,p;double cond,detA=1;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin(0.2*i)-0.64*exp(0.1/i);b=0.16;c=-0.064;YaSuoJZ(C,a,b, c); //将矩阵A中元素压缩存储在C中λm1=mifa(C); //对A用幂法求出模最大的特征值λm1λs=fanmifa(C); //对A用反幂法求出模最小的特征值λsYaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中for(j=0;j<=500;j++) //对A进行平移，平移量为λm1，平移后矩阵元素压缩存储在C中C[2][j]=C[2][j]-λ?m1;λm=mifa(C);λm2=λm1+λm; //λm1与λm2是矩阵的最大最小特征值if(λm1>λm2) //判断A最大最小特征值{λ501=λm1;λ1=λm2;}else{λ501=λm2;λ1=λm1;}printf("数值分析计算实习第一题\n\n ZY1515105 樊雪松\n\n （1）A的最大最小以及模最小的特征值\n");printf("A的最小特征值λ1=%.13e\n",λ1);printf("A的最大特征值λ501=%.13e\n",λ501);printf("A的模最小特征值λs=%.13e\n",λs);printf("\n（2）与数μk最接近的特征值\n");printf("\t要求接近的值\t\t\t实际求得的特征值\n");YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中for(k=1;k<=39;k++){p=λ1+k*(λ501-λ1)/40;for(j=0;j<=501;j++)C[2][j]=C[2][j]-p;λ=fanmifa(C)+p;printf("μ%d=%.13e λ%d=%.13e\n",k,p,k,λ);YaSuoJZ(C,a,b, c); //还原矩阵A中元素并压缩存储在C中}printf("\n（3）计算A的条件数cond(A)和行列式detA\n");cond=λm1/λs;daizhuangLU(C);for(j=0;j<=500;j++)detA=detA*C[2][j];printf("A的条件数cond(A)=%.13e\n",cond);printf("A的行列式detA=%.13e\n",detA);getch();}三、运行结果数值分析计算实习第一题ZY1515105 樊雪松（1）A的最大最小以及模最小的特征值A的最小特征值λ1=-1.0700113615018e+001A的最大特征值λ501=9.7246340987773e+000A的模最小特征值λs=-5.5579107942295e-003（2）与数μk最接近的特征值要求接近的值实际求得的特征值μ1=-1.0189494922173e+001 λ1=-1.0182934033146e+001 μ2=-9.6788762293280e+000 λ2=-9.5857074250676e+000 μ3=-9.1682575364831e+000 λ3=-9.1726724239280e+000 μ4=-8.6576388436383e+000 λ4=-8.6522840078976e+000 μ5=-8.1470201507934e+000 λ5=-8.0934838086753e+000 μ6=-7.6364014579485e+000 λ6=-7.6594054076924e+000 μ7=-7.1257827651036e+000 λ7=-7.1196846486912e+000 μ8=-6.6151640722588e+000 λ8=-6.6117643393973e+000 μ9=-6.1045453794139e+000 λ9=-6.0661032265951e+000 μ10=-5.5939266865690e+000 λ10=-5.5851010526284e+000 μ11=-5.0833079937241e+000 λ11=-5.1140835298122e+000 μ12=-4.5726893008792e+000 λ12=-4.5788721768651e+000 μ13=-4.0620706080344e+000 λ13=-4.0964709262599e+000 μ14=-3.5514519151895e+000 λ14=-3.5542112157508e+000 μ15=-3.0408332223446e+000 λ15=-3.0410900181333e+000 μ16=-2.5302145294997e+000 λ16=-2.5339703111304e+000 μ17=-2.0195958366549e+000 λ17=-2.0032307695635e+000μ18=-1.5089771438100e+000 λ18=-1.5035576112274e+000μ19=-9.9835845096511e-001 λ19=-9.9355860600754e-001μ20=-4.8773975812023e-001 λ20=-4.8704267388496e-001μ21=2.2878934724645e-002 λ21=2.2317362495748e -002μ22=5.3349762756952e-001 λ22=5.3241747420686e -001μ23=1.0441163204144e+000 λ23=1.0528989626935e+000μ24=1.5547350132593e+000 λ24=1.5894458818809e+000μ25=2.0653537061042e+000 λ25=2.0603304602743e+000μ26=2.5759723989490e+000 λ26=2.5580755970728e+000μ27=3.0865910917939e+000 λ27=3.0802405093071e+000μ28=3.5972097846388e+000 λ28=3.6136208676923e+000μ29=4.1078284774837e+000 λ29=4.0913785104506e+000μ30=4.6184471703285e+000 λ30=4.6030353782791e+000μ31=5.1290658631734e+000 λ31=5.1329242838984e+000μ32=5.6396845560183e+000 λ32=5.5949063480833e+000μ33=6.1503032488632e+000 λ33=6.0809338570269e+000μ34=6.6609219417080e+000 λ34=6.6803540921116e+000μ35=7.1715406345529e+000 λ35=7.2938774481266e+000μ36=7.6821593273978e+000 λ36=7.7171117142356e+000μ37=8.1927780202427e+000 λ37=8.2252200140502e+000μ38=8.7033967130876e+000 λ38=8.6486660651935e+000μ39=9.2140154059324e+000 λ39=9.2542003445750e+000（3）计算A 的条件数cond(A)和行列式detAA 的条件数cond(A)=1.9252042739022e+003A 的行列式detA=2.7727861417521e+118四、结果分析设A 的n 个线性无关的特征向量为1x ，2x ,…,n x ，其相对应的特征值满足的关系为n λλλλ≥≥≥> 321。

【2018最新】数值分析考试试卷-精选word文档 (43页)

本文部分内容来自网络整理，本司不为其真实性负责，如有异议或侵权请及时联系，本司将立即删除！== 本文为word格式，下载后可方便编辑和修改！ ==数值分析考试试卷篇一：数值分析试题及答案一、单项选择题（每小题3分，共15分）1. 3.142和3.141分别作为?的近似数具有（）和（）位有效数字. A．4和3 B．3和2 C．3和4 D．4和4?2dx?12. 已知求积公式1f?x?6f?1??Af(213)?6f(2)，则A＝（）1112A． 6 B．3C．2 D．33. 通过点?x0,y0?,?x1,y1?的拉格朗日插值基函数l0?x?,l1?x?满足（）A．＝0，l1?x1??0B．l0?x0?＝0，l1?x1??1C．l0?x0?＝1，l1?x1??1D．l0?x0?＝1，l1?x1??14. 设求方程f?x??0的根的牛顿法收敛，则它具有（）敛速。

A．超线性 B．平方 C．线性 D．三次?x?1?2x2?x3?0?2x1?2x2?3x3?35. 用列主元消元法解线性方程组???x1?3x2 ?2作第一次消元后得到的第3个方程（ A．B．?2x2?1.5x3?3.5C．?2x2?x3?3D．x2?0.5x3??1.5单项选择题答案1.A2.D3.D4.C5.B二、填空题（每小题3分，共15分）.）1. 设X?(2,3,?4), 则||X||1?，||X||2?2. 一阶均差f?x0,x1??TC0x?3?3. 已知n?3时，科茨系数4. 因为方程内有根。

?18,C1?3??C2?3??38，那么C3?3??f?x??0f?x??x?4?2?0在区间?1,2?上满足，所以在区间5. 取步长h?0.1，用欧拉法解初值问题y??y??y?2x??y?1??1?的计算公式 .填空题答案1三、计算题（每题15分，共60分）y?1?x的一组数据：21. 已知函数求分f?1.5?。

数值分析第五版计算实习题

数值分析计算实习题第二章2-1程序：clear;clc;x1=[0.2 0.4 0.6 0.8 1.0];y1=[0.98 0.92 0.81 0.64 0.38];n=length(y1);c=y1(:);for j=2:n %求差商for i=n:-1:jc(i)=(c(i)-c(i-1))/(x1(i)-x1(i-j+1));endendsyms x df d;df(1)=1;d(1)=y1(1);for i=2:n %求牛顿差值多项式df(i)=df(i-1)*(x-x1(i-1));d(i)=c(i)*df(i);enddisp('4次牛顿插值多项式');P4=vpa(collect((sum(d))),5) %P4即为4次牛顿插值多项式,并保留小数点后5位数pp=csape(x1,y1, 'variational');%调用三次样条函数q=pp.coefs;disp('三次样条函数');for i=1:4S=q(i,:)*[(x-x1(i))^3;(x-x1(i))^2;(x-x1(i));1];S=vpa(collect(S),5)endx2=0.2:0.08:1.08;dot=[1 2 11 12];figureezplot(P4,[0.2,1.08]);hold ony2=fnval(pp,x2);x=x2(dot);y3=eval(P4);y4=fnval(pp,x2(dot));plot(x2,y2,'r',x2(dot),y3,'b*',x2(dot),y4,'co');title('4次牛顿插值及三次样条');结果如下：4次牛顿插值多项式P4 = - 0.52083*x^4 + 0.83333*x^3 - 1.1042*x^2 + 0.19167*x + 0.98三次样条函数x∈[0.2,0.4]时， S = - 1.3393*x^3 + 0.80357*x^2 - 0.40714*x + 1.04 x∈[0.4,0.6]时，S = 0.44643*x^3 - 1.3393*x^2 + 0.45*x + 0.92571x∈[0.6,0.8]时，S = - 1.6964*x^3 + 2.5179*x^2 - 1.8643*x + 1.3886 x∈[0.8,1.0]时，S =2.5893*x^3 - 7.7679*x^2 + 6.3643*x - 0.80571输出图如下2-3（1）程序：clear;clc;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];%插值点n=length(y1);a=ones(n,2);a(:,2)=-x1';c=1;for i=1:nc=conv(c,a(i,:));endq=zeros(n,n);r=zeros(n,n+1);for i=1:n[q(i,:),r(i,:)]=deconv(c,a(i,:));%wn+1/(x-xk) end 三次样条插值曲线4次牛顿插值曲线Dw=zeros(1,n);for i=1:nDw(i)=y1(i)/polyval(q(i,:),x1(i));%系数endp=Dw*q;syms x L8;for i=1:nL8(i)=p(n-i+1)*x^(i-1);enddisp('8次拉格朗日插值');L8=vpa(collect((sum(L8))),5)xi=0:64;yi=polyval(p,xi);figureplot(xi,yi,x1,y1,'r*');hold ontitle('8次拉格朗日插值');结果如下：8次拉格朗日插值L8 =- 3.2806e-10*x^8 + 6.7127e-8*x^7 - 5.4292e-6*x^6 + 0.00022297*x^5 - 0.0049807*x^4 + 0.060429*x^3 - 0.38141*x^2 + 1.3257*x输出图如下：第五章4-1（3）程序：clc;clear;y= @(x) sqrt(x).*log(x);a=0;b=1;tol=1e-4;p=quad(y,a,b,tol);fprintf('采用自适应辛普森积分结果为： %d \n', p);结果如下：采用自适应辛普森积分结果为： -4.439756e-01第九章9-1(a)程序：clc;clear;a=1;b=2;%定义域h=0.05;%步长n=(b-a)/h;y0=1;%初值f= @(x,y) 1/x^2-y/x;%微分函数Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值%以下根据改进欧拉公式求解for i=1:nxn=Xn(i);xnn=Xn(i+1);yn=Yn(i);yp=yn+h*f(xn,yn);yc=yn+h*f(xnn,yp);yn=(yp+yc)/2;Yn(i+1)=yn;endXn=Yn;%以下根据经典四阶R-K法公式求解for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;enddisp(' 改进欧拉法四阶经典R-K法'); disp([Xn' Yn'])结果如下：改进欧拉法四阶经典R-K法1 10.99887 0.998850.99577 0.99780.99114 0.996940.98532 0.996340.97857 0.996030.97111 0.996060.96311 0.996450.9547 0.997230.94598 0.998410.93705 10.92798 1.0020.91883 1.00440.90964 1.00730.90045 1.01060.89129 1.01430.88218 1.01840.87315 1.02290.86421 1.02780.85538 1.03310.84665 1.0388（b）程序：clc;clear;a=0;b=1;%定义域H=[0.1 0.025 0.01];%步长y0=1/3;%初值f= @(x,y) -50*y+50*x^2+2*x;%微分函数xi=linspace(a,b,11);Y=1/3*exp(-50*xi)+xi.^2;%准确解Ym=zeros(1,11);for j=1:3h=H(j);n=(b-a)/h;Xn=linspace(a,b,n+1);%将定义域分为n等份Yn=zeros(1,n);%结果矩阵Yn(1)=y0;%赋初值for i=1:nxn=Xn(i);yn=Yn(i);k1=f(xn,yn);k2=f(xn+h/2,yn+h/2*k1);k3=f(xn+h/2,yn+h/2*k2);k4=f(xn+h,yn+h*k3);yn=yn+h/6*(k1+2*k2+2*k3+k4);Yn(i+1)=yn;endfor k=1:11m=0.1/h;Ym(k)=Yn(1+(k-1)*m);enddelta=Ym-Y;fprintf('步长为： %d \n', h);disp(' 四阶经典R-K法准确解误差'); disp([Ym' Y' delta'])end结果如下：步长为： 1.000000e-01四阶经典R-K法准确解误差0.33333 0.33333 04.6055 0.012246 4.593263.062 0.040015 63.022864.05 0.09 863.9611844 0.16 118431.6235e+05 0.25 1.6235e+052.2256e+06 0.36 2.2256e+063.0509e+07 0.49 3.0509e+074.1823e+08 0.64 4.1823e+085.7333e+09 0.81 5.7333e+097.8594e+10 1 7.8594e+10步长为： 2.500000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.013015 0.012246 0.000768940.040063 0.040015 4.82e-050.090037 0.09 3.6857e-050.16004 0.16 3.6723e-050.25004 0.25 3.6722e-050.36004 0.36 3.6722e-050.49004 0.49 3.6722e-050.64004 0.64 3.6722e-050.81004 0.81 3.6722e-051 1 3.6722e-05步长为： 1.000000e-02四阶经典R-K法准确解误差0.33333 0.33333 00.012256 0.012246 9.5673e-060.040016 0.040015 7.8252e-070.090001 0.09 6.6347e-070.16 0.16 6.6226e-070.25 0.25 6.6225e-070.36 0.36 6.6225e-070.49 0.49 6.6225e-070.64 0.64 6.6225e-070.81 0.81 6.6225e-071 1 6.6225e-07由结果可知，步长越小，结果越精确。

数值分析计算实习题（二）

数值分析计算实习题（二）数值分析计算实习题（二）SY1004114 全昌彪一：算法设计方案概述：本题采用fortran90语言编写程序，依据题目要求，采用带双步位移QR分解法求出所给矩阵的所有特征值，并求出相应于其实特征值的特征向量，以及相关需要给出的中间结果。

1、矩阵的A的初始化（赋值）：利用子函数initial（a,n）来实现，返回n×n 维二维数组a。

2、A矩阵的拟上三角化：利用子函数hessenberg(a,n)，在对矩阵进行QR分解前进行拟上三角化，这样可以提高计算效率，减少计算量，返回A矩阵的相似矩阵Hessenberg阵A(n-1)。

3、对A(n-1)进行带双步位移QR分解得出Ｃｍ及A矩阵的所有特征值，这一步利用了两个子函数eigenvalue(a,n,lamda,lamdai)和qrresolve(b,c,m)带双步位移QR分解可以加速收敛。

每次QR分解前先进行判断，若可以直接得到矩阵的特征值，则对矩阵直接降阶处理；若不可以，则进行QR分解，这样就进一步减少了计算量，提高了计算效率。

考虑到矩阵A可能有复特征值，采用两个一维数组lamda(n)及lamdai(n)分别存储其实部和虚部。

在双步位移处理及降阶过程中，被分解的矩阵Ａk(m ×m)及中间矩阵M k(m×m)的维数随m不断减少而降阶，于是引入了动态矩阵C(m×m)和B(m×m)分别存储，在使用前，先声明分配内存，使用结束后立即释放内存。

返回A(n-1)经双步位移QR分解后的矩阵及Ａ矩阵的所有特征值。

4、特征向量的求解：采用子函数eigenvector(a,lamda)实现求解A矩阵的属于实特征值的特征向量。

核心算法为高斯列主元消去法，(A-λI)x=b,b=0,回代过程令x(10)=1，即可求出对应于每一实特征值的特征向量的各个元素。

5、相关输出结果：所有数据均采用e型输出，数据保留到小数点后12位。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

矩阵与数值分析2018秋上机实习
1. 用秦九韶算法编程计算f (x )=1+x +x 2+⋯+x 50在x =1.00001处的值。

2. 设f (x )=54x 6+45x 5−102x 4−69x 3+35x 2+16x −4.在区间[−2,2] 上画出函数，
（1）使用Newton 迭代法找出该区间上的5个根，并计算e i+1/e i 2和e i+1/e i ，由此判断哪个根是1阶收敛，哪个根是2阶收敛？（2）使用割线法计算这5个根，并判断哪个根是线性收敛，哪个是超线性收敛？
3. 令H 表示n ×n 的Hilbert 矩阵，其中(i,j)元素是1/(i +j −1), b 是元素全为1的向量，用Gauss 消去法求解Hx =b,其中取(a) n =2; (b) n =5; (c) n =10.
4. 已知方程组
[ 3−1−13−1⋱⋱⋱ −13−1−13]
[x 1⋮x n ]=[ 21⋮12] 分别用Jacobi 迭代和Gauss-Seidel 迭代求解方程组，精确到小数点后6位
5. 用共轭梯度法求解第3题中的方程组
6. 令f (x )=e |x|,x ∈[−1,1]，分别用等距节点和Chebyshev 的零点去插值f(x)，等距节点包括左右两个端点，分别取n =5,10,15,20，画出插值函数以及原函数的图并比较，观察有没有龙格现象发生。

7. 编程求解教材183页例3，并计算出样条函数在插值结点及相邻结点的中点处的导数值，并画出原函数及插值函数，原函数的导函数及插值函数的导函数的图像。

把步长变为0.1重复上述操作。

8. （1）给定数据点(x i ,x i 2),x i =0,1n ,2n ,…,1，当n =5,10,15,20,25,30时分别用直线拟合这组数据点并注意观察当点数逐渐增加时直线的表达式的变化.（2）计算函数f (c 1,c 2)=∫(x 2−c 1−c 2x )2dx 1
0的最小值，并解释与（1）的关系
9. 已知常微分方程
{du dx =2x u +x 2e x x ∈[1,2],u (1)=0
，分别用Euler 法，改进的Euler 法，Runge-Kutta 法去求解该方程，步长选为0.1,0.05,0.01.画图观察求解效果。

要求：
1.考试前提交作业（word 形式提交，包括代码和实验结果），主题写“学号+姓名+学部(学院)”发送至邮箱zhuke_2015@ ，文件名“学号+姓名+学部学院”。

考试后提交的数值试验部分成绩记为零分。

2.可用任何一种语言编程。