数值分析计算实习题(二)
数值分析计算实习题二
《数值分析》计算实习题二算法设计方案1．主要计算步骤：计算函数f(x,y)在拟合所需的节点处的函数值。
将各拟合节点（x i，y j）分别带入非线性方程组0.5 cos t + u + v + w – x = 2.67t + 0.5 sin u + v + w – y = 1.070.5t + u + cos v + w – x =3.74t + 0.5u + v + sin w – y =0.79解非线性方程组得解向量（t ij，u ij，v ij，w ij）。
对数表z(t,u)进行分片二次代数插值，求得对应（t ij，u ij）处的值，即为f(x i，y j) 的值。
对上述拟合节点分别进行x，y最高次数为k（k=0,1,2,3…）次的多项式拟合。
每次拟合后验证误差大小，直到满足要求。
2．求解非线性方程组选择Newton迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle分解法。
3．对z(t，u)进行插值选择分片二次插值。
4．拟合基函数φr(x)ψs(y)选择为φr(x)=x r，ψs(y)=y s。
拟合系数矩阵c通过连续两次解线性方程组求得。
一．源程序#include "stdio.h"#include "stdlib.h"#include "math.h"void Doolittle(double *A,int n,int *M)//功能说明：对n阶矩阵A进行选主元的Doolittle分解//参数说明：A：欲进行分解的方阵，同时也是返回参数，分解后的结果// 存储于A中// n：方阵A的维数// M；（返回参数）n维向量，记录选主元过程中行交换的次序{int i,j,k,t;double *s;double Maxs,temp;s=(double*) calloc(n,sizeof(double));for(k=0;k<n;k++){for(i=k;i<n;i++){s[i]=A[i*n+k];for(t=0;t<k;t++) s[i]-= A[i*n+t] * A[t*n+k];}Maxs=abs(s[k]); M[k]=k;for(i=k+1;i<n;i++){if(Maxs<abs(s[i])){Maxs=abs(s[i]);M[k]=i;}}if(M[k]!=k){for(t=0;t<n;t++){temp=A[k*n+t];A[k*n+t]=A[M[k]*n+t];A[M[k]*n+t]=temp;}temp=s[k];s[k]=s[M[k]];s[M[k]]=temp;}if(Maxs<(1e-14)){s[k]=1e-14;printf("%.16e方阵奇异\n",Maxs);}A[k*n+k]=s[k];for(j=k+1;(j<n)&&(k<n-1);j++){for(t=0;t<k;t++) A[k*n+j]-=A[k*n+t]*A[t*n+j];A[j*n+k]=s[j]/A[k*n+k];}}}void Solve_LUEquation(double* A,int n,double* b,double* x) //功能说明：解方程LUx＝b，其中L、U共同存储在A中//参数说明：A：经Doolittle分解后的方阵// n：方阵A的维数// b：方程组的右端向量// x：（返回参数）方程组的解向量{int i,t;for(i=0;i<n;i++){x[i]=b[i];for(t=0;t<i;t++) x[i]-=A[i*n+t]*x[t];}for(i=n-1;i>-1;i--){for(t=i+1;t<n;t++) x[i]-=A[i*n+t]*x[t];x[i]/=A[i*n+i];}}void Transpose(double *A,int m,int n,double* AT)//功能说明：求m×n阶矩阵A的转置AT//参数说明：A：已知m×n阶矩阵// m：A的行数// n：A的列数// AT：（返回参数）A的转置矩阵（n×m）{int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++) AT[j*m+i]=A[i*n+j];}void Solve_LEquation(double* A,int n,double* B,double* x,int m) //功能说明：解线性方程组Ax＝B，该函数可对系数矩阵相同// 而右端向量不同的多个方程组同时求解。
数值分析实验报告2
实验报告实验项目名称函数逼近与快速傅里叶变换实验室数学实验室所属课程名称数值逼近实验类型算法设计实验日期班级学号姓名成绩512*x^10 - 1280*x^8 + 1120*x^6 - 400*x^4 + 50*x^2 - 1并得到Figure，图像如下：实验二：编写程序实现[-1,1]上n阶勒让德多项式，并作画（n=0,1，…，10 在一个figure中）。
要求：输入Legendre(-1,1,n)，输出如a n x n+a n-1x n-1+…多项式。
在MATLAB的Editor中建立一个M-文件，输入程序代码，实现勒让德多项式的程序代码如下：function Pn=Legendre(n,x)syms x;if n==0Pn=1;else if n==1Pn=x;else Pn=expand((2*n-1)*x*Legendre(n-1)-(n-1)*Legendre(n-2))/(n);endx=[-1:0.1:1];A=sym2poly(Pn);yn=polyval(A,x);plot (x,yn,'-o');hold onend在command Windows中输入命令：Legendre（10），得出的结果为：Legendre(10)ans =(46189*x^10)/256 - (109395*x^8)/256 + (45045*x^6)/128 - (15015*x^4)/128 + (3465*x^2)/256 - 63/256并得到Figure，图像如下：实验三：利用切比雪夫零点做拉格朗日插值，并与以前拉格朗日插值结果比较。
在MATLAB的Editor中建立一个M-文件，输入程序代码，实现拉格朗日插值多项式的程序代码如下：function [C,D]=lagr1(X,Y)n=length(X);D=zeros(n,n);D(:,1)=Y';for j=2:nfor k=j:nD(k,j)=(D(k,j-1)- D(k-1,j-1))/(X(k)-X(k-j+1));endendC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(X(k)));m=length(C);C(m)= C(m)+D(k,k);end在command Windows 中输入如下命令：clear,clf,hold on;k=0:10;X=cos(((21-2*k)*pi)./22); %这是切比雪夫的零点Y=1./(1+25*X.^2);[C,D]=lagr1(X,Y);x=-1:0.01:1;y=polyval(C,x);plot(x,y,X,Y,'.');grid on;xp=-1:0.01:1;z=1./(1+25*xp.^2);plot(xp,z,'r')得到Figure ，图像如下所示：比较后发现，使用切比雪夫零点做拉格朗日插值不会发生龙格现象。
北航数值分析计算实习题目二 矩阵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 λλ---++=的两个根。
计算方法与实习答案1-2
绪论
习题1——10：设 f ( x) = 8 x 5 − 0.4 x 4 + 4 x 3 − 9 x + 1 用秦九韶法求f(3)。 解：
8 − 0.4
24 8 23.6
0
−9
1
x=3
70.8 74.8
224.4 224.4
673.2 664.2
1992.6 1993.6
∴ f(3)=1993.6
第一章 绪论 练习
1．《计算方法》课程主要研究以计算 机为工具的 数值 分析方法 ，并评价 该算法的计算误差。 2．近似值作四则运算后的绝对误差限 公式为 ε ( x1 − x2 ) ≤ ε ( x1 ) + ε ( x2 ) ，近似值 1.0341的相对误差限不大于 1 ×10−2 ， 则它至少有三位有效数字。 4
ln(103 ) ∴k ≥ ln(2) ≥ 9.965
2 2 2
∴需二分10次 需二分 次
方程求根——二分法
习题2——2：用二分法求方程2e-x-sinx=0在区 间[0,1]内的1个实根，要求3位有效数字。
解：1）判断是否在该区间有且仅有一个根 f(0)=2>0，f(1)=2/e-sin1≈-0.1<0， f’(x)=-2e-x-cosx，f’=-3,-2/e-cos1<0 2）判断二分次数 由(b-a)/2k+1=1/2k+1≤1/2*10-3，解得k≥3ln10/ln2≥9.965， 所以需要二分10次，才能满足精度要求。
∴ x≈2.981
方程求根
f (xk )(xk − xk −1) xk +1 = xk − f (xk ) − f (xk −1)
习题2——11：用割线法求方程x3-2x-5=0的根，要 求精确到4位有效数字，取x0=2, x1=2.2。
数值分析计算实习题
数值分析计算实习题-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)%观察图像，类似抛物线，故用二次多项式拟合。
数值分析上机实习题及答案.docx
方詡文金兴:爭［数值分析］2017-2018第二学期上机实习题1：编程计算亍丄，其中C= 4. 4942x10307,给出并观察计算结心C"果，若有问题，分析之。
解:mat lab 编程如下:E) funct ion diy i ti formatlong g;n 二input C 输入ii 值= c= 4.4942E307; sum 二0; s 二 0;E3 for i = l:n s = l/ (c*i);>> diyiti 输入n 值:10 104.6356e-308 >> diyiti输入ri 值:1001004.6356e-308 >> diyiti 输入n 值:1000 10004.6356e-308 >> diyiti揄入n 值* 1000001000004・ 6356e-308 >> diyiti输入n 值；1000000001000000004.6356e-308图二：运行结果Mat lab 中，forma t long g 对双精度，显示15位定点或浮点格式，由上图 可知，当输入较小的n 值5分别取10, 100, 1000, 100000, 100000000）的时候, 结果后面的指数中总是含有e-308,这和题目中的C 值很相似，我认为是由于分 母中的C 值相对于n 值过大，出现了 “大数吃小数”的现彖，这是不符合算法原 则的。
2：利用牛顿法求方程-1^ = 2于区间241的根，考虑不同初值下牛顿法的收敛情况。
解：牛顿法公式为:利用mat lab 编程function di2ti21 3i=l ；2 2.85208156699784 xO 二input （'输入初值x0：‘ ）;A 二[i x0];3 2.55030468822809 t=x0+ (x0-log (xO) -2) /(1-1/xO) ; %迭代函数4 1.91547247100476 三 while (abs (t _x0)>0.01)i=i+l; 5 0.37867158538991 xO 二 t; 6 0.774964959780275 A = [A;i xO];t =x0+(x0-log(xO)-2)/(1-1/xO): 7 4.11574081641933 cnd| 8 5.04162436446126 disp (A);96.81782826645596当输入初值二3的时候并不能收敛。
数值分析计算实习题
插值法1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64] 上作图.(1)用这9 个点作8 次多项式插值Ls(x).(2)用三次样条( 第一边界条件)程序求S(x).从得到结果看在[0,64] 上，哪个插值更精确；在区间[0,1] 上，两种插值哪个更精确解：(1) 拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64]; y1=[0 1 2 3 4 5 6 7 8]; n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:n l=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l; endLs=simplify(Ls) % 为所求插值多项式Ls(x).输出结果为Ls = -/*xA2+95549/72072*x-1/00*xA8-2168879/0*xA4+19/0*xA7+657859/*xA3+33983/ 0*xA5-13003/00*xA6(2) 三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=s plin e(x1,y1,x2);p=po Iyfit(x2,y3,3); % 得到三次样条拟合函数S=p(1)+p(2)*x+ p(3)*x^2+p(4)*xA3 % 得到S(x) 输出结果为：S =/6464-2399/88*x+/1984*xA2+2656867/624*xA3⑶ 在区间［0,64］上，分别对这两种插值和标准函数作图,Plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y="X函数曲线，红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64] 上三次样条插值更精确。
数值分析第二次作业答案
练习1 已知410=x，211=x，432=x。
（1）推导以这3点作求积节点在[0，1]上的插值求积公式；（2）指明该求积公式所具有的代数精度； （3）用所求的公式计算dxx ⎰12解：按题设原式是插值型的，故有32434121414321100=⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-=⎰dx x x A同样，容易计算出3202==A A ，于是有求积公式⎪⎭⎫⎝⎛+⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛≈⎰433221314132)(1f f f dx x f由于原式含有3个节点，按定理1它至少有2阶精度。
考虑到其对称性，可以猜到它可能有3阶精度。
事实上，对于3)(x x f =原式左右两端相等。
此外，容易验证原式对4)(x x f =不准确，故所构造的求积公式确实有3阶精度。
（3）31]43221412[31222102=⎪⎭⎫⎝⎛⨯+⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛⨯⨯≈⎰dx x31432141214341101-=⎪⎭⎫ ⎝⎛-⎪⎭⎫⎝⎛-⎪⎭⎫ ⎝⎛-⎪⎭⎫ ⎝⎛-=⎰dx x x A2. 取7个等距节点（包括区间端点）分别用复化梯形公式和复化辛甫生公式求积分2lnxdx的近似值（取6位小数）解：（1）复化梯形公式])(2)()([2)(11∑⎰-=++=≈n k k ban x f b f a f h T dx x f385139.0])(2)2()1([1211616=++=∴∑-=k k x f f f T(2)复化辛甫生公式])(2)(4)()([6)(11021∑∑⎰-=-=++++=≈n k k n k k n bax f xf b f a f h S dx x f ∴ ])(2)(4)2()1([3161212213∑∑==++++⨯=k k k k x f xf f f S≈0.386 287而 38629436.0ln21=⎰xdx3. 用梯形格式求解初值问题⎩⎨⎧=≤<++-='2)1(6.,y x x y y )1(1 1 ，（取步长h ＝0.2，小数点后至少保留6位） 解：梯形格式为)],(),([2111+++++=n n n n n n y x f y x f h y y ，于是⇒++-+++-+=+++ 1 1 ,)]()[(2111n n n n n n x y x y h y y),(222112 +++++-=++n n n n x x hh y hh y,2,1,0=n取步长h ＝0.2，由初值20=y 计算得147709.2)6.1(069422.2)4.1(018182.2)2.1(321=≈=≈=≈y y y y y y4. 对初值问题⎩⎨⎧=>=+'1)0(00y x y y , 试证明用欧拉预-校格式所求得的近似解为,2,1,022, )-(1=+=n hh y nn (其中h 为步长)证明： ,2,1,0)],(),([2),(1111 =⎪⎩⎪⎨⎧++=+=++++n y x f y x f hy y y x hf y y n n n n n n n n n n 将y y x f -=) ( ,代入，于是有⎪⎩⎪⎨⎧--+=-=+++)(2)1(111n n n n n n y y hy y y h y 整理后，有)-(1n n y hh y 221+=+反复递推得 )-(101212y hh y n n +++=由1)0(0==y y ，故得,2,1,022, )-(1=+=n hh y nn。
西南交通大学数值分析上机实习
目录解题： (1)题目一： (1)1.1计算结果 (1)1.2结果分析 (1)题目二： (2)2.1计算结果 (2)2.2结果分析 (3)题目三： (4)3.1计算结果 (4)3.2结果分析 (5)总结 (5)附录 (6)Matlab程序： (6)题目一： (6)第一问Newton法： (6)第二问Newton法： (6)第一问Steffensen加速法： (7)第二问Steffensen加速法： (7)题目二 (8)1、Jacobi迭代法 (8)2、Causs-Seidel迭代法 (8)题目三： (9)题目一：分别用牛顿法，及基于牛顿算法下的Steffensen 加速法(1)求ln(x +sin x )=0的根。
初值x0分别取0.1, 1,1.5, 2, 4进行计算。
(2)求sin x =0的根。
初值x0分别取1，1.4，1.6, 1.8，3进行计算。
分析其中遇到的现象与问题。
1.1计算结果求ln(x +sin x )=0的根，可变行为求解x-sinx-1=0的根。
1.2结果分析从结果对比我们可发现牛顿—Steffensen 加速法比牛顿法要收敛的快，牛顿法对于初值的选取特别重要，比如第（1）问中的初值为4的情况，100次内没有迭代出来收敛解，而用Steffensen 加速法，7次迭代可得；在第（2）问中的初值为1.6的情况，收敛解得31.4159，分析其原因应该是x x f cos )('=，x0=1.62π≈,0)('≈x f ；迭代式在迭代过程中会出现分母趋近于0，程序自动停止迭代的情况，此时得到的x 往往非常大，而在第一问中我们如果转化为用x+sinx=1,则可以收敛到结果。
用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。
(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T,(2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,b2=[5,0,-10]T,(3)A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T2.1计算结果初值均为0矩阵带入(1)A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T2) A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,b2=[5,0,-10]TT2.2结果分析ρ小于1,故方程组雅可比迭代收第一小题的经计算谱半径为5427B(=).0敛。
数值分析参考答案(第二章)doc资料
证明：
（1）
得证。
+
得证。
14． 求 及 。
解：
若
则
15．证明两点三次埃尔米特插值余项是
解：
若 ，且插值多项式满足条件
插值余项为
由插值条件可知
且
可写成
其中 是关于 的待定函数，
现把 看成 上的一个固定点，作函数
根据余项性质，有
由罗尔定理可知，存在 和 ，使
即 在 上有四个互异零点。
根据罗尔定理， 在 的两个零点间至少有一个零点，
数值分析参考答案(第二章)
第二章插值法
1．当 时， ,求 的二次插值多项式。
解：
则二次拉格朗日插值多项式为
2．给出 的数值表
X
0.4
0.5
0.6
0.7
0.8
lnx
-0.916291
-0.693147
-0.510826
-0.356675
-0.223144
用线性插值及二次插值计算 的近似值。
解：由表格知，
若采用线性插值法计算 即 ，
则
若采用二次插值法计算 时，
3．给全 的函数表，步长 若函数表具有5位有效数字，研究用线性插值求 近似值时的总误差界。
解：求解 近似值时，误差可以分为两个部分，一方面，x是近似值，具有5位有效数字，在此后的计算过程中产生一定的误差传播；另一方面，利用插值法求函数 的近似值时，采用的线性插值法插值余项不为0，也会有一定的误差。因此，总误差界的计算应综合以上两方面的因素。
解：函数 的 展式为
其中
又 是次数为 的多项式
为 阶多项式
为 阶多项式
依此过程递推，得 是 次多项式
数值分析计算实习二
void atra(double a[N][N], int n);//矩阵转置
void aadd(double a[N], double b[N], double c[N], int n, int co);//向量加法
void aaadd(double a[N][N], double b[N][N], double c[N][N], int n, int co);//矩阵加法
结果：
1、 拟上三角矩阵 如下。
2、迭代最终矩阵如下，迭代次数为14次。
3、矩阵所有特征值如下。
从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。
4、对应于实特征值的特征向量如下
问题发现与讨论
1、双步位移QR迭代过程中迭代矩阵 始终为拟上三角矩阵。
2、双步位移QR方法与基本QR方法比较
amul(wr, ur, D, m+1);
aaadd(A, D, A, m+1, -1);
amul(ur, pr, D, m+1);
aaadd(A, D, A, m+ห้องสมุดไป่ตู้, -1);
}
}
}
k++;
if (k > 10000)//算法失败判定
{
printf("算法失败\n");
break;
}
}
for (int i = 0; i < N; i++)
doubleaneiji(double a[N], double b[N], int n);//向量内积
void aamul(double a[N], double b[N][N], double c[N], int n);//矩阵向量乘法
数值分析题库及标准答案
模 拟 试 卷（一）一、填空题（每小题3分，共30分）1．有3个不同节点的高斯求积公式的代数精度是 次的.2．设152210142-⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦A ，342⎛⎫ ⎪=- ⎪ ⎪⎝⎭x ，则 ∞A = .， 1x = ______.3．已知y =f (x )的均差（差商）01214[,,]3f x x x =，12315[,,] 3f x x x =，23491[,,]15f x x x =，0238[,,] 3f x x x =, 那么均差423[,,]f x x x = .4．已知n =4时Newton －Cotes 求积公式的系数分别是：,152,4516,907)4(2)4(1)4(0===C C C 则)4(3C ＝ .5．解初始值问题00(,)()y f x y y x y '=⎧⎨=⎩的改进的Euler 方法是 阶方法；6．求解线性代数方程组123123123530.13260.722 3.51x x x x x x x x x --=⎧⎪-++=⎨⎪++=⎩的高斯—塞德尔迭代公式为 ，若取(0)(1,1,1)=-r x, 则(1)=rx .7．求方程()x f x =根的牛顿迭代格式是 .8．01(), (),, ()n x x x l l L l 是以整数点01, ,, ,n x x x L 为节点的Lagrange 插值基函数，则()nk jk k x x =∑l= .9．解方程组=Ax b 的简单迭代格式(1)()k k +=+xBx g 收敛的充要条件是 .10．设(-1)1,(0)0,(1)1,(2)5f f f f ====，则()f x 的三次牛顿插值多项式为 ，其误差估计式为 .二、综合题（每题10分，共60分）1．求一次数不超过4次的多项式()p x 满足：(1)15p =，(1)20p '=，(1)30p ''=(2)57p =，(2)72p '=.2．构造代数精度最高的形式为10101()()(1)2xf x dx A f A f ≈+⎰的求积公式，并求出其代数精度.3．用Newton 法求方程2ln =-x x 在区间) ,2(∞内的根, 要求8110--<-kk k x x x .4．用最小二乘法求形如2y a bx =+的经验公式拟合以下数据：5．用矩阵的直接三角分解法解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡71735 30103421101002014321x x x x .6 试用数值积分法建立求解初值问题0(,)(0)y f x y y y '=⎧⎨=⎩的如下数值求解公式1111(4)3n n n n n hy y f f f +-+-=+++，其中(,),1,,1i i i f f x y i n n n ==-+.三、证明题（10分）设对任意的x ，函数()f x 的导数()f x '都存在且0()m f x M '<≤≤,对于满足20Mλ<<的任意λ，迭代格式1()k k k x x f x λ+=-均收敛于()0f x =的根*x .参考答案一、填空题1．5； 2. 8, 9 ; 3.9115； 4. 1645； 5. 二； 6. (1)()()123(1)(1)()213(1)(1)(1)312(330.1)/5(220.7)/6(12)*2/7k k k k k k k k k x x x x x x x x x ++++++⎧=++⎪=+-⎨⎪=--⎩, (0.02，0.22，0.1543) 7. 1()1()k k k k k x f x x x f x +-=-'-; 8. j x ; 9. ()1B ρ<;10.32(4)11,()(1)(1)(2)/24(1,2)66x x x f x x x x ξξ+-+--∈-二、综合题 1．差商表：233234()1520(1)15(1)7(1)(1)(2)5432p x x x x x x x x x x =+-+-+-+--=++++其他方法：设233()1520(1)15(1)7(1)(1)()p x x x x x ax b =+-+-+-+-+ 令(2)57p =，(2)72p '=，求出a 和b. 2．取()1,f x x =，令公式准确成立，得：0112A A +=, 011123A A +=, 013A =, 116A =. 2()f x x =时，公式左右14=；3()f x x =时，公式左15=, 公式右524=∴ 公式的代数精度2=. 3．此方程在区间) ,2(∞内只有一个根s ，而且在区间（2，4）内。
数值分析计算实习作业二
数值分析计算实习题二学号：姓名：院系：2015年11月28日一、算法采用拟上三角化的双步位移的QR分解法计算特征值，并利用列主元素Guass消去法求实数特征值对应的特征向量，定义全局变量n=10。
1 主程序在定义一维和二维数组时将数组的大小定为n+1（大小为11），这样在使用单个元素时的角标可以与矩阵或向量元素的角标相对应。
先输入矩阵A的各个元素；保留数组A，定义数组B的元素与A相同，即B=A，后面的矩阵变换用数组B来进行；对B进行拟上三角化；输出拟上三角化后的矩阵B；对拟上三角化后的B进行QR分解，输出矩阵Q、R、RQ；对拟上三角化的数组矩阵B进行QR分解；输出A的全部特征值；判断所有特征值，如果特征值为实数（即储存特征值虚部的数组中元素为0），利用列主元素高斯消去法求其特征向量；输出实特征值对应的特征向量。
2 拟上三角化程序对于r=1,2，……，n-2循环执行（1）求第r列的第r+2个元素到第n个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r的循环。
如果和不为0，则继续执行。
（2）计算第r列从第r+1到第n个元素的平方和的平方根dr=当B[r+1][r]小于等于0时，取cr=dr，否则取cr=-drhr=cr*cr-cr*B[r+1][r]（3）令(0,...,0,[1][],[2][],...,[][])Tur B r r cr B r r B n r=+-+（4）计算()/()()/()()()/()()()()()()()T T Tpr B ur hr qr B ur hr tr pr ur hr wr qr tr ur B B wr ur ur pr ====-=--（5）执行下一次循环。
3 拟上三角化后的矩阵的QR 分解初始令Q=I 。
对于r=1,2，……，n-1循环执行（1）求第r 列的第r+1个元素到第n 个元素的平方和，如果和为0，则跳出该次循环直接进行下一次r 的循环。
数值实验暑期实习题目
数值实验2011暑期实习题目********************************************************************实验一：迭代法求方程根的收敛性实验目的：体会在非线性方程求根的迭代法中，迭代函数和初始迭代值的选取对迭代收敛性的影响。
实验内容：考虑一个简单的代数方程324100xx +-=，针对上述方程，可以构造多种迭代法，如下列几种迭代格式：①321410n n n x x x x +=--++；②13211(10)2n n x x +=-；③12110()4n x x +=+；④321241038n n n n n nx x x x x x ++-=-+。
实验要求：（1）取定某个初始值0x ，按如上四种迭代格式进行计算，它们的收敛性如何？重复选取不同的初始值，反复实验，并自行设计一种比较形象的记录方式（如利用Matlab 的图形功能等），分析四种迭代法的收敛性与初值选取的关系。
（2）选取第④种迭代格式（Newton 迭代法），取不同的初始值进行迭代，结果如何？并分析迭代法对不同的初值是否有差异。
（3）对上述四种迭代格式，编制Steffensen 迭代程序，选取不同的初始值，输出迭代次数和方程的根，并与（1）、（2）中的结果进行比较。
********************************************************************实验二：多项式插值的振荡现象，即插值的龙格现象实验目的：体会在用多项式插值方法逼近一个函数时，使用的节点越多，插值多项式的次数就越高，随着插值多项式次数的增加，插值函数是否也更加逼近给定的函数。
实验内容：考虑龙格函数21()125f x x=+在区间[1,1]-的一个等距划分，分点取为 21,0,1,2,...,i ix i n n=-+=则拉格朗日插值多项式为201()()125nn ii iL x l x x ==+∑其中的()(0,1,2,,)i l x i n =是n 次拉格朗日插值基函数。
数值分析上机实习题
数值分析上机实习题第2章插值法1. 已知函数在下列各点的值为试⽤四次⽜顿插值多项式）（x p 4及三次样条韩式)(S x (⾃然边界条件)对数据进⾏插值。
⽤图给出(){}10,11,1,0,08.02.0,,x i =+=i x y i i ，）（x p 4及)(x S Python 代码import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.font_manager import FontPropertiesfont_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc",size=12) #求⽜顿n 次均差 def qiujuncha(x,f,n): for i in range(1,n): for j in range(4,i-1,-1):f[j]= (f[j] - f[j-1])/(x[j]-x[j-i]) #根据⽜顿多项式求值 def niudun(x,f,x1): sum = f[0]; tmp = 1;for i in range(1,5): tmp *= (x1-x[i-1]) sum = sum + f[i]*tmp return sum#⽜顿插值画图 def drawPic(x,f):x1 = np.linspace(0.2, 1, 100) plt.plot(x1, niudun(x,f,x1))plt.title(u"⽜顿四次插值",fontproperties=font_set) plt.xlabel(u"x 轴",fontproperties=font_set) plt.ylabel(u"y 轴", fontproperties=font_set) plt.show() def qiu_h(x,h): n = len(x) -1 for i in range(n): print(i)h[i] = x[i+1]-x[i]#⾃然边界条件下的三次样条插值求Mdef qiu_m(h,f,o,u,d):n = len(h)o[0] = 0u[n] = 0d[n] = d[0] = 0a = []for i in range(1,n):u[i] = h[i-1]/(h[i-1]+h[i])for i in range(1,n):o[i] = h[i]/(h[i-1]+h[i])for i in range(1,n-1):d[i] = 6*(f[i+1]-f[i])/(h[i-1]+h[i])t = [0 for i in range(5)]t[0] =2t[1] = o[0]a.append(t)for i in range(1,n):t = [0 for i in range(5)]t[i - 1] = u [i + 1]t[i] = 2t[i + 1] = o [i + 1]a.append(t)t = [0 for i in range(5)]t[n - 1] = u[n]t[n] = 2a.append(t)tmp = np.linalg.solve(np.array(a),np.array(d))m = []for i in range(5):m.append(tmp[i])return m#根据三次条插值函数求值def yangtiao(x1,m,x,y,h,j):returnm[j]*(x[j+1]-x1)**3/(6*h[j])+m[j+1]*(x1-x[j])**3/(6*h[j])+(y[j]-m[j]*h[j]**2/6)*(x[j+1]-x1)/h[j] +(y[j+1]-m[j+1]*h[j]**2/6)*(x1-x[j])/h[j] def main():x = [0.2, 0.4, 0.6, 0.8, 1.0]y = [0.98, 0.92, 0.81, 0.64, 0.38]f = y[:]f1 = y[:]h = [0.2,0.2,0.2,0.2]u = [0 for n in range(5)]d = [0 for n in range(5)]o = [0 for n in range(5)] qiujuncha(x,f,4) qiujuncha(x,f1,2)m = qiu_m(h,f1,o,u,d) x1 = np.linspace(0.2, 0.4, 10)p1= plt.plot(x1, yangtiao(x1,m,x,y,h,0),color='red') x1 = np.linspace(0.4, 0.6, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 1),color='red') x1 = np.linspace(0.6, 0.8, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 2),color='red') x1 = np.linspace(0.8, 1.0, 10)plt.plot(x1, yangtiao(x1, m, x, y, h, 3),color='red') x1 = np.linspace(0.2, 1.0, 40)p2 = plt.plot(x1,niudun(x,f,x1),color='green') plt.xlabel(u"x 轴", fontproperties=font_set) plt.ylabel(u"y 轴",fontproperties=font_set) plt.title("三次样条插值和⽜顿插值")plt.legend(labels=[u'三次样条插值',u'⽜顿插值'],prop=font_set,loc="best") plt.show() main()实验结果运⾏结果可得插值函数图（如图1-1），4次⽜顿插值函数）（x p 4和三次样条插值函数)(x S 如下：)6.0(*)4.0(*)2.0(625.0)4.0(*)2.0(*3.098.0)(4-------=x x x x x x x P 98.0)8.0(*)6.0(*)4.0(*)2.0(*20833.0+-----x x x x]4.0,2.0[),2.0(467.4)4.0(9.4)2.0(167.1)(S 3∈-+-+-=x x x x x]6.0,4.0[),4.0(113.4)6.0(6467.4)4.0(575.1)6.0(167.1)(S 33∈-+-+----=x x x x x x ]8.0,6.0[),6.0(2.3)8.0(113.4)6.0(575.1)(S 3∈-+-+--=x x x x x]0.1,8.0[),8.0(9.1)0.1(2.3)(S ∈-+-=x x x x图1-1三次样条插值和⽜顿插值图2.在区间[-1,1]上分别取n = 10,20⽤两组等距节点对龙格函数做多项式插值三次样条插值，对每个n值画出插值函数及图形。
数值分析(第五版)计算实习题第四章作业
第四章：1、（1）：复合梯形建立m文件：function t=natrapz(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b); f=feval(fname,a+h:h:b-h+0.001*h); t=h*(0.5*(fa+fb)+sum(f));输入：>> syms x>> f=inline('sqrt(x).*log(x);'); >> natrapz(f,eps,1,10)输出：ans =-0.417062831779470输入：>> syms x>> f=inline('sqrt(x).*log(x);'); >> natrapz(f,eps,1,100)输出：ans =-0.443117908008157输入：>> syms x>> f=inline('sqrt(x).*log(x);'); >> natrapz(f,eps,1,1000)输出：ans =-0.444387538997162复合辛普森建立m文件：function t=comsimpson(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0.001*h);f2=feval(fname,a+h/2:h:b-h+0.001*h);t=h/6*(fa+fb+2*sum(f1)+4*sum(f2));输入：>> syms x>> f=inline('sqrt(x).*log(x);');>> format long；>>comsimpson(f,eps,1,10)输出：ans =-0.435297890074689输入：>>syms x>>f=inline('sqrt(x).*log(x);');>>comsimpson(f,eps,1,100)输出：ans =-0.444161178415673输入：>>syms x>>f=inline('sqrt(x).*log(x);');>>comsimpson(f,eps,1,1000)输出：ans =-0.444434117614180（2）龙贝格建立m文件：function [RT,R,wugu,h]=Romberg(fun,a,b,wucha,m) %RT是龙贝格积分表%R是数值积分值%wugu是误差估计%h是最小步长%fun是被积函数%a b是积分下、上限%m是龙贝格积分表中行最大数目%wucha是两次相邻迭代值的绝对误差限n=1;h=b-a;wugu=1;x=a;k=0;RT=zeros(4,4);RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2;while((wugu>wucha)&(k<m)|(k<4))k=k+1;h=h/2;s=0;for j=1:nx=a+h*(2*j-1);s=s+feval(fun,x);endRT(k+1,1)=RT(k,1)/2+h*s;n=2*n;for i=1:kRT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);endwugu=abs(RT(k+1,k)-RT(k+1,k+1));endR=RT(k+1,k+1);输入：>>fun=inline('sqrt(x).*log(x)');>> [RT,R,wugu,h]=Romberg(fun,eps,1,1e-5,13)输出：RT =1 至5 列-0.000000268546145 0 0 0-0.245064670140209 -0.326752804004897 0 0-0.358104125949240 -0.395783944552250 -0.400386020588741 0 0-0.408090073087781 -0.424752055467295 -0.426683262861631 -0.427100679405645 0-0.429474601629505 -0.436602777810080 -0.437392825966266 -0.437562819031419 -0.437603847029951-0.438389494461832 -0.441361125405941 -0.441678348578999 -0.441746372747455 -0.4417627788404596 列-0.441766844267449R =-0.441766844267449wugu =4.065426989774412e-06h =0.031250000000000（3）自适应辛普森输入：>> f=inline('sqrt(x).*log(x)');>> q=quad(f,0,1,1e-4)输出：q =-0.4439755729517282.（1）复合辛普森建立m文件function q=combinesimpson2(F,x0,a,b,n)%复合Simpson多元求积公式%F—被积函数%x0—被积函数自变量%[a,b]积分区间%n—区间份数x=linspace(a,b,n+1);q=0;for k=1:nq=q+subs(F,x0,x(k))+4*subs(F,x0,(x(k)+x(k+1))/2)+subs(F,x0,x(k+1)); endq=q*(b-a)/n/6;输入：>> clear>> syms x y;>> F=exp(-x.*y);>> s=combinesimpson2(combinesimpson2(F,'x',0,1,4),'y',0,1,4)输出：s =exp(-1)/576 + exp(-1/2)/144 + exp(-1/4)/72 + exp(-3/4)/144 + exp(-1/8)/36 +exp(-3/8)/36 + exp(-5/8)/72 + exp(-7/8)/72 + (5*exp(-1/16))/144 + exp(-3/16)/24 + exp(-5/16)/36 + exp(-7/16)/36 + exp(-9/16)/144 + exp(-1/32)/36 + exp(-3/32)/18 + exp(-5/32)/36 + exp(-7/32)/36 + exp(-9/32)/36 + exp(-15/32)/36 + exp(-21/32)/36 + exp(-1/64)/36 + exp(-3/64)/18 + exp(-5/64)/18 + exp(-7/64)/18 + exp(-9/64)/36 + exp(-15/64)/18 + exp(-21/64)/18 + exp(-25/64)/36 + exp(-35/64)/18 + exp(-49/64)/36 + 47/576>> double(s)ans =0.796599967946203高斯求积公式function q=gaussquad(F,x0,a,b,n)%Gauss求积公式%F—被积函数%x0—被积函数自变量%[a,b]积分区间%n—节点个数syms t;F=subs(F,x0,(b-a)/2*t+(a+b)/2);[x,A]=gausspoints(n);q=(b-a)/2*sum(A.*subs(F,t,x));输入：>> clear>> syms x y;F=exp(-x.*y);>> s=gaussquad(gaussquad(F,x,0,1,4),y,0,1,4)输出：s =0.7966（2）复合辛普森输入：>> syms x y;>> f=exp(-x.*y);>> s=combinesimpson2(combinesimpson2(f,y,0,sqrt(1-x^2),4),x,0,1,4)输出：s =(3^(1/2)*(exp(-3^(1/2)/4) + 2*exp(-3^(1/2)/8) + 2*exp(-3^(1/2)/16) + 2*exp(-(3*3^(1/2))/16) + 4*exp(-3^(1/2)/32) + 4*exp(-(3*3^(1/2))/32) + 4*exp(-(5*3^(1/2))/32) + 4*exp(-(7*3^(1/2))/32) + 1))/576 + (7^(1/2)*(exp(-(3*7^(1/2))/16) + 2*exp(-(3*7^(1/2))/32) + 2*exp(-(3*7^(1/2))/64) + 2*exp(-(9*7^(1/2))/64) + 4*exp(-(3*7^(1/2))/128) + 4*exp(-(9*7^(1/2))/128) + 4*exp(-(15*7^(1/2))/128) + 4*exp(-(21*7^(1/2))/128) + 1))/1152 + (15^(1/2)*(exp(-15^(1/2)/16) + 2*exp(-15^(1/2)/32) + 2*exp(-15^(1/2)/64) + 2*exp(-(3*15^(1/2))/64) + 4*exp(-15^(1/2)/128) + 4*exp(-(3*15^(1/2))/128) + 4*exp(-(5*15^(1/2))/128) + 4*exp(-(7*15^(1/2))/128) + 1))/1152 + (15^(1/2)*(exp(-(7*15^(1/2))/64) + 2*exp(-(7*15^(1/2))/128) + 2*exp(-(7*15^(1/2))/256) + 2*exp(-(21*15^(1/2))/256) + 4*exp(-(7*15^(1/2))/512) + 4*exp(-(21*15^(1/2))/512) + 4*exp(-(35*15^(1/2))/512) + 4*exp(-(49*15^(1/2))/512) + 1))/1152 + (39^(1/2)*(exp(-(5*39^(1/2))/64) + 2*exp(-(5*39^(1/2))/128) + 2*exp(-(5*39^(1/2))/256) + 2*exp(-(15*39^(1/2))/256) + 4*exp(-(5*39^(1/2))/512) + 4*exp(-(15*39^(1/2))/512) + 4*exp(-(25*39^(1/2))/512) + 4*exp(-(35*39^(1/2))/512) + 1))/1152 + (55^(1/2)*(exp(-(3*55^(1/2))/64) + 2*exp(-(3*55^(1/2))/128) + 2*exp(-(3*55^(1/2))/256) + 2*exp(-(9*55^(1/2))/256) + 4*exp(-(3*55^(1/2))/512) + 4*exp(-(9*55^(1/2))/512) + 4*exp(-(15*55^(1/2))/512) + 4*exp(-(21*55^(1/2))/512) + 1))/1152 + (63^(1/2)*(exp(-63^(1/2)/64) + 2*exp(-63^(1/2)/128) + 2*exp(-63^(1/2)/256) + 2*exp(-(3*63^(1/2))/256) + 4*exp(-63^(1/2)/512) + 4*exp(-(3*63^(1/2))/512) + 4*exp(-(5*63^(1/2))/512) + 4*exp(-(7*63^(1/2))/512) + 1))/1152 + 1/24>> double(s)ans =0.670113633359095。
数值分析计算实习题
1.下列数据点的插值x 0 1 4 9 16 25 36 49 64y 0 1 2 3 4 5 6 7 8可以得到平方根函数的近似，在区间[0,64]上作图.(1)用这9个点作8次多项式插值Ls(x).(2)用三次样条(第一边界条件)程序求S(x).从得到结果看在[0,64]上，哪个插值更精确；在区间[0,1]上，两种插值哪个更精确？解：(1)拉格朗日插值多项式，求解程序如下syms x l;x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];n=length(x1);Ls=sym(0);for i=1:nl=sym(y1(i));for k=1:i-1l=l*(x-x1(k))/(x1(i)-x1(k));endfor k=i+1:nl=l*(x-x1(k))/(x1(i)-x1(k));endLs=Ls+l;endLs=simplify(Ls) %为所求插值多项式Ls(x).输出结果为Ls =-/*x^2+95549/72072*x-1/00*x^8-2168879/0*x^4+19/0*x^7+657859/*x^3+33983/ 0*x^5-13003/00*x^6(2)三次样条插值，程序如下x1=[0 1 4 9 16 25 36 49 64];y1=[0 1 2 3 4 5 6 7 8];x2=[0:1:64];y3=spline(x1,y1,x2);p=polyfit(x2,y3,3); %得到三次样条拟合函数S=p(1)+p(2)*x+p(3)*x^2+p(4)*x^3 %得到S(x)输出结果为：S =/6464-2399/88*x+/1984*x^2+2656867/624*x^3(3)在区间[0,64]上，分别对这两种插值和标准函数作图，plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')蓝色曲线为y=函数曲线,红色曲线为拉格朗日插值函数曲线，黄色曲线为三次样条插值曲线可以看到蓝色曲线与黄色曲线几乎重合，因此在区间[0,64]上三次样条插值更精确。
数值分析计算实习题（二）
数值分析计算实习题（二）数值分析计算实习题（二）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位。
数值分析实习第二题
for(j = 0; j < Max_N; j++)
{
cache[i] += Matrix[j][i]*Vect[j];
}
}
for(i = 0; i < Max_N; i++)
{
Vect_Ret[i] = cache[i];
}
}
void VECT_EVA(char Max_N, double Vect_new[], double Vect_bef[])
{
u[k] = 0;
}
}
VECT_MULTI(N, u, w, 1/h);
MAT_T_VEC(N, A, w, p);
MAT_VEC(N, A, w, q);
VECT_MULTI(N, u, w, -VECT_T_VEC(N, p, w));
VECT_PLUS(N, q, w, w);
for(k = 0 ;k < N; k++)
//向量转置乘向量
{
char n;
double power_norm;
for(n = 0, power_norm = 0; n < Max_N; n++)
{
power_norm += Vector_T[n]*Vector[n];
}
return(power_norm);
}
void VECT_PLUS(char Max_N, double Vector_1[], double Vector_2[], double Vector_Ret[])
MAT_VEC(Max_N, A, w, q);
VECT_MULTI(Max_N, u, w, -VECT_T_VEC(Max_N, p, w));
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析计算实习题（二）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位。
对于Ａ矩阵的拟上三角化结果集双步位移QR分解结果比较庞大，为了更好的显示，还采用了ｆ型输出，保留５位小数。
所有计算精度水平Ｅ＝10-12。
二：算法fortran源程序!此函数用于给A赋值，即初始化A矩阵subroutine initial(a,n)integer :: i,jdimension a(n,n)double precision ado i=1,ndo j=1,nif (i==j) thena(i,j)=1.5*cos(i+1.2*j)elsea(i,j)=sin(0.5*i+0.2*j)end ifend doend doend subroutine initial!此函数用于对Ａ矩阵进行拟上三角化subroutine hessenberg(a,n)implicit noneinteger::i,j,r,ndimension a(n,n),u(n),p(n),q(n),w(n) double precision a,u,p,q,w,d,t,c,h,z,e e=1.0d-12 !精度控制outer1:do r=1,n-2inner1:do i=1,np(i)=0.0q(i)=0.0w(i)=0.0end do inner1d=0.0t=0.0z=0inner2:do i=r+2,nif (a(i,r)/=0) thenz=1end ifend do inner2if (z==1) theninner3:do i=r+1,nd=sqrt(d**2+a(i,r)**2)end do inner3if (a(r+1,r)>0) thenc=-delsec=dend ifh=c**2-c*a(r+1,r)inner4:do i=1,nif (i<=r) thenu(i)=0else if (i==r+1) thenu(i)=a(r+1,r)-celseu(i)=a(i,r)end ifend do inner4inner5:do i=1,ndo j=1,np(i)=p(i)+u(j)*a(j,i)/hq(i)=q(i)+a(i,j)*u(j)/hend doend do inner5inner6:do i=1,nt=t+p(i)*u(i)end do inner6t=t/hinner7:do i=1,nw(i)=q(i)-t*u(i)end do inner7inner8:do i=1,ndo j=1,na(i,j)=a(i,j)-w(i)*u(j)-u(i)*p(j)end doend do inner8end ifend do outer1do i=1,ndo j=1,nif (abs(a(i,j))<=e) then!a(i,j)输出精度控制a(i,j)=0end ifend doend dowrite(*,*) ('矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)（位有效数字e型输出）:') write(*,100) ((a(i,j),j=1,n),i=1,n)write(*,*) ('矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)(保留五位小数f型输出):') write(*,200) ((a(i,j),j=1,n),i=1,n)!输出格式控制100format (10e20.12)200format (10f9.5)end subroutine hessenberg!QR分解子函数，求出Ａ（ｎ－１）矩阵QR分解的结果!主要是用于eigenvalue(a,n,lamda,lamdai)的调用subroutine qrresolve(b,c,m)implicit noneinteger::i,j,r,mdimension b(m,m),c(m,m),u(m),p(m),q(m),w(m),v(m)double precision b,c,u,p,q,w,v,d,c1,h,t,zouter1:do r=1,m-1inner1:do i=1,m!中间量初始化p(i)=0.0q(i)=0.0w(i)=0.0v(i)=0.0end do inner1d=0.0t=0.0z=0inner2:do i=r+1,mif (b(i,r)/=0) thenz=1end ifend do inner2if (z==1) theninner3:do i=r,md=sqrt(d**2+b(i,r)**2) end do inner3if (b(r,r)>0) thenc1=-delsec1=dend ifh=c1**2-c1*b(r,r)inner4:do i=1,mif (i<r) thenu(i)=0else if (i==r) thenu(i)=b(r,r)-c1elseu(i)=b(i,r)end ifend do inner4inner5:do i=1,mdo j=1,mv(i)=v(i)+u(j)*b(j,i)/h end doend do inner5inner6:do i=1,mdo j=1,mb(i,j)=b(i,j)-u(i)*v(j) end doend do inner6inner7:do i=1,mdo j=1,mp(i)=p(i)+u(j)*c(j,i)/hq(i)=q(i)+c(i,j)*u(j)/hend doend do inner7inner8:do i=1,mt=t+p(i)*u(i)end do inner8t=t/hinner9:do i=1,mw(i)=q(i)-t*u(i)end do inner9inner10:do i=1,mdo j=1,mc(i,j)=c(i,j)-w(i)*u(j)-u(i)*p(j)end doend do inner10end ifend do outer1end subroutine qrresolve!此函数用于求A矩阵的特征值（包含实特征值即复特征值）!核心算法为双步位移QR方法（中间调用了QR分解法子函数）!此程序最大的特点是只用了一个goto语句!（按课本上给出个goto语句将更易编写，笔者也编写了运行结果一样）!为了清晰显示，每一个循环语句，if语句基本上都有命名subroutine eigenvalue(a,n,lamda,lamdai)implicit nonedimension a(n,n),lamda(n),lamdai(n)double precision,allocatable,dimension(:,:) :: b,c !动态数组，用于保存Mk及C矩阵元素，维数是变化的integer :: i,j,k,m,n,Ldouble precision a,s,t,temp1,temp2,temp3,deta,lamda,lamdai,EE=1.0D-12 !精度m=nL=150 !循环迭代最大次数限定lamda=0lamdai=0do k=1,Louter1:if (abs(a(m,m-1))<=E) thenlamda(m)=a(m,m)lamdai(m)=0m=m-1 !m自减，降一维10 inner1:if (m==1) thenlamda(m)=a(1,1)lamdai(m)=0exitelse if (m==0) thenexitelsecycleend if inner1else!Dk子块求根temp1=(a(m-1,m-1)+a(m,m))/2.0temp2=a(m-1,m-1)*a(m,m)-a(m,m-1)*a(m-1,m)deta=temp1**2-temp2inner2:if (deta>0) thenlamda(m)=temp1+sqrt(deta)lamdai(m)=0lamda(m-1)=temp1-sqrt(deta)lamdai(m-1)=0else if (deta<0) thenlamda(m)=temp1lamda(m-1)=temp1lamdai(m)=sqrt(-deta)lamdai(m-1)=-sqrt(-deta)elselamda(m)=temp1lamda(m-1)=temp1lamdai(m)=0lamdai(m-1)=0end if inner2inner3: if (m==2) thenexitelseinner4: if (abs(a(m-1,m-2))<=E) thenm=m-2 !m自减，降两维goto 10elseinner5: if (k==L) thenexitelseallocate(b(m,m)) !动态数组分配内存声明allocate(c(m,m))inner6: do i=1,mdo j=1,mc(i,j)=a(i,j)end doend do inner6s=a(m-1,m-1)+a(m,m)t=a(m-1,m-1)*a(m,m)-a(m,m-1)*a(m-1,m) b=matmul(c,c) !fotran90可直接用矩阵乘法inner7: do i=1,mdo j=1,mb(i,j)=b(i,j)-s*c(i,j)if (i==j) thenb(i,j)=b(i,j)+tend ifend doend do inner7call qrresolve(b,c,m) !调用qr分解子函数inner8: do i=1,mdo j=1,ma(i,j)=c(i,j)end doend do inner8 !结束，释放可分配内存deallocate(b)deallocate(c)end if inner5end if inner4end if inner3end if outer1end dodo i=1,ndo j=1,nif (abs(a(i,j))<=e) then!a(i,j)输出精度控制a(i,j)=0end ifend doend dowrite(*,*) ('矩阵A(n-1)双步位移QR分解的结果（位有效数字e型输出）:') write(*,100) ((a(i,j),j=1,n),i=1,n)write(*,*) ('矩阵A(n-1)双步位移QR分解的结果（保留五位小数f型输出）:') write(*,200) ((a(i,j),j=1,n),i=1,n)write(*,*) ('矩阵A的全部特征值:')do i=1,nwrite(*,300) lamda(i),'+',lamdai(i),'i'end do!输出格式控制语句100format (10e20.12)200format (10f9.5)300format (1x,1e20.12,a2,1x,1e20.12,a1)end subroutine eigenvalue!此函数用于求矩阵A属于实特征值的特征向量!核心算法为高斯列主元消去法，(A-λI)x=b,b=0,x(10)=1回代subroutine eigenvector(a,lamda)parameter n=10dimension a(n,n),m(n),b(n),x(n),e(n,n),a_max(n),i_k(n) double precision a,lamda,b,x,t,a_max,sum,e,i_kinteger m,i,j,k,lcall initial(a,n) !调用子函数，对A进行初始化!单位矩阵do i=1,ne(i,i)=1end do!各参数初始化b=0m=0a_max=0t=0a=a-lamda*e !A-λI!消元过程do k=1,n-1do i=k,nif (abs(a(i,k))>=a_max(k)) thena_max(k)=a(i,k) !选列主元m(k)=iend ifend dodo j=k,nt=a(k,j)a(k,j)=a(m(k),j)a(m(k),j)=tt=b(k)b(k)=b(m(k))b(m(k))=tend dodo i=k+1,ni_k(i)=a(i,k)/a(k,k)do j=k+1,na(i,j)=a(i,j)-i_k(i)*a(k,j)end dob(i)=b(i)-i_k(i)*b(k)end doend do!回代过程x(n)=1do k=n-1,1,-1sum=0do j=k+1,nsum=sum+a(k,j)*x(j)end dox(k)=(b(k)-sum)/a(k,k)end dowrite (*,100) 'A的属于lamda=',lamda,'的特征向量:' write (*,200) x(:)100format (3x,a14,1e19.12,a11)200format (1x,1e20.12)end subroutine eigenvector!主函数，用于调用各个子函数program main_functionimplicit nonedimension a(10,10),lamda(10),lamdai(10),x(10)integer :: i,j,ndouble precision a,lamda,lamdai,xn=10call initial(a,n) !A矩阵初始化call hessenberg(a,n) !A(n-1)矩阵拟上三角化call eigenvalue(a,n,lamda,lamdai) !求A矩阵的所有特征值!求A矩阵属于实特征值的特征向量call eigenvector(a,lamda(1))call eigenvector(a,lamda(4))call eigenvector(a,lamda(5))call eigenvector(a,lamda(8))call eigenvector(a,lamda(9))call eigenvector(a,lamda(10))end program main_function三：程序输出结果汇总：矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)（12位有效数字e型输出）:-0.882751733065E+00 -0.993313316411E-01 -0.110334955595E+01 -0.760044676830E+00 0.154909623840E+00-0.194659168629E+01 -0.878226438751E-01 -0.925593370941E+00 0.603254375819E+00 0.151882370093E+00-0.234787833029E+01 0.237236973917E+01 0.181929100012E+01 0.323780721082E+00 0.220580321938E+000.210269272044E+01 0.181611739616E+00 0.127884456394E+01 -0.638050158135E+00 -0.415402766732E+000.000000000000E+00 0.172827457015E+01 -0.117146801290E+01 -0.124383914708E+01 -0.639976054127E+00-0.200283235265E+01 0.292496128495E+00 -0.641286862951E+00 0.978337580632E-01 0.255770647387E+000.000000000000E+00 0.000000000000E+00 -0.129166965220E+01 -0.111160388379E+01 0.117134648592E+01-0.130735637674E+01 0.180370952121E+00 -0.424641232157E+00 0.798856288534E-01 0.160878411848E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.156012616609E+01 0.812504527179E+000.442174836453E+00 -0.358869742717E-01 0.469176590186E+00 -0.273656195621E+00 -0.735907113693E-010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 -0.770777431331E+00-0.158305209469E+01 -0.304282234567E+00 0.252871461915E+00 -0.670993649410E+00 0.254459252711E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00-0.746342926455E+00 -0.270872578791E-01 -0.948652518297E+00 0.119586353771E+00 0.192915386004E-010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 -0.770179463681E+00 -0.469763787838E+00 0.498822338663E+00 0.113768224333E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.701315606266E+00 0.158220658011E+00 0.386262135075E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.484380788945E+00 0.399280281752E+00矩阵A拟上三角化得到的Hessenberg矩阵A(n-1)(保留五位小数f型输出):-0.88275 -0.09933 -1.10335 -0.76004 0.15491 -1.94659 -0.08782 -0.92559 0.60325 0.15188-2.34788 2.37237 1.81929 0.32378 0.22058 2.10269 0.18161 1.27884 -0.63805 -0.415400.00000 1.72827 -1.17147 -1.24384 -0.63998 -2.00283 0.29250 -0.64129 0.09783 0.255770.00000 0.00000 -1.29167 -1.11160 1.17135 -1.30736 0.18037 -0.42464 0.07989 0.160880.00000 0.00000 0.00000 1.56013 0.81250 0.44217 -0.03589 0.46918 -0.27366 -0.073590.00000 0.00000 0.00000 0.00000 -0.77078 -1.58305 -0.30428 0.25287 -0.67099 0.254460.00000 0.00000 0.00000 0.00000 0.00000 -0.74634 -0.02709 -0.94865 0.11959 0.019290.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -0.77018 -0.46976 0.49882 0.113770.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.70132 0.15822 0.386260.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.48438 0.39928矩阵A(n-1)双步位移QR分解的结果（12位有效数字e型输出）:0.338303924721E+01 0.894875525475E+00 -0.895678244508E+00 0.846512632866E-01 0.261267793959E+000.161043313340E+01 -0.102255900702E+01 0.937223229328E-01 -0.100257738993E+01 -0.408625114186E+000.000000000000E+00 -0.211848218734E+01 -0.236152888518E+01 -0.345573222905E-01 -0.473663572352E-010.181640642861E+01 -0.231836979306E+00 -0.143548828808E+00 -0.653703238456E+00 0.322737505939E-010.000000000000E+00 0.355512159597E+00 -0.252851058727E+01 -0.637526428692E+00 0.202388052178E-010.183862769119E+01 0.186937662765E+00 -0.293254867876E+00 0.198706730246E+01 0.100463175221E+010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.157754818551E+01 -0.139621667385E-01-0.697199695782E+00 0.155545542881E+00 0.840521541031E-02 -0.815438441765E-01 -0.108612039515E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 -0.148403977622E+01-0.100530588729E+00 0.424954671617E-01 0.262346756601E-01 0.104004221001E+00 -0.118074774235E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00-0.694860436802E+00 -0.528901458582E+00 0.267916098808E+00 -0.596236908619E+00 -0.491158600319E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+000.178846387945E+00 -0.126620169756E+01 0.471455617468E-01 0.290478168004E+00 -0.357041410381E-010.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.935588710088E+00 0.187740659272E+00 0.136024981245E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.636050766399E+00 0.273703606786E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.000000000000E+000.000000000000E+00 0.000000000000E+00 0.000000000000E+00 0.245568085231E-04 0.565162119254E-01矩阵A(n-1)双步位移QR分解的结果（保留五位小数f型输出）:3.38304 0.89488 -0.89568 0.08465 0.26127 1.61043 -1.02256 0.09372 -1.00258 -0.408630.00000 -2.11848 -2.36153 -0.03456 -0.04737 1.81641 -0.23184 -0.14355 -0.65370 0.032270.00000 0.35551 -2.52851 -0.63753 0.02024 1.83863 0.18694 -0.29325 1.987071.004630.00000 0.00000 0.00000 1.57755 -0.01396 -0.69720 0.15555 0.00841 -0.08154 -0.108610.00000 0.00000 0.00000 0.00000 -1.48404 -0.10053 0.04250 0.02623 0.10400 -0.118070.00000 0.00000 0.00000 0.00000 0.00000 -0.69486 -0.52890 0.26792 -0.59624 -0.491160.00000 0.00000 0.00000 0.00000 0.00000 0.17885 -1.26620 0.04715 0.29048 -0.035700.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.93559 0.18774 0.136020.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.63605 0.273700.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00002 0.05652矩阵A的全部特征值:0.338303924721E+01 + 0.000000000000E+00i-0.232349638731E+01 + -0.893040543164E+00i-0.232349638731E+01 + 0.893040543164E+00i0.157754818551E+01 + 0.000000000000E+00i-0.148403977622E+01 + 0.000000000000E+00i-0.980531067180E+00 + -0.113949139468E+00i-0.980531067180E+00 + 0.113949139468E+00i0.935588710088E+00 + 0.000000000000E+00i0.565046144245E-01 + 0.000000000000E+00i0.636062363900E+00 + 0.000000000000E+00iA的属于lamda= 0.338303924721E+01的特征向量: -0.436739834602E+00-0.906357312937E+00-0.196334018248E+01-0.108285729157E+01-0.126929785310E+01-0.107244858764E+010.362510824633E+000.168290589916E+010.211274357722E+010.100000000000E+01A的属于lamda= 0.157754818551E+01的特征向量: 0.621734685009E+00-0.111511734070E+00-0.248377525794E+01-0.130686124886E+01-0.381560716897E+010.811730907979E+01-0.123917001440E+01-0.680029414194E+000.269189633346E+010.100000000000E+01A的属于lamda=-0.148403977622E+01的特征向量: -0.173077810726E+020.240818399765E+020.413384989796E+00-0.857202156181E+010.928742697832E-01-0.783266447227E-01-0.637425984951E+00-0.340318975687E+00-0.378486114544E+000.100000000000E+01A的属于lamda= 0.935588710088E+00的特征向量: 0.279242914534E+010.159824395408E+01-0.520736954810E+00-0.166788993965E+01-0.122571220310E+020.724123969831E+01-0.539823439712E+010.284100862311E+02-0.121651397115E+020.100000000000E+01A的属于lamda= 0.565046144245E-01的特征向量:-0.510500561103E+01-0.488632284214E+010.950516391275E+01-0.678833013260E+00-0.960433257896E+01-0.304574395607E+010.157487434454E+02-0.739503192243E+01-0.710966608732E+010.100000000000E+01A的属于lamda= 0.636062363900E+00的特征向量:0.474503250282E+010.315786893686E+010.172995116830E+02-0.198004944728E+01-0.318752699418E+020.779400967078E+01-0.100425554520E+020.167075308672E+020.131053073601E+020.100000000000E+01请按任意键继续. . .四：编写程序过程中发现的一个问题本程序运行正确，但是在编写过程中曾经产生一个很小的错误，将矩阵M k=A k2-s A k+t I 表达式中错写成了M k=A k2-s A k-t I，依旧可以得到Ａ矩阵的所有特征值及属于实特征值的特征向量，而且和本程序的答案完全吻合，只是Ａ(n-1)进行带双步位移QR分解结束后得到的矩阵C m相同（因为和同学对比了计算结果，发现不对），于是我调试了很久（也重新写过带双步位移QR分解子程序）方才发现这个小变动，而这个小变动只影响带双步位移QR分解，并不影响特征值的求解。