北航数值分析报告大作业第八题
北航 数值分析第二次大作业(带双步位移的QR方法)
一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。
总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。
为了防止矩阵在后面的计算中被破坏保存A[][]。
2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。
由于矩阵的QR 分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。
这里用函数void QuasiTriangularization()来实现，函数形参为double型N维方阵double a[][N]。
3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。
用子程序void QR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。
4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。
为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。
若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。
这里用函数intTwoStepDisplacement_QR()来实现。
这是用来存储计算得到的特征值的二维数组。
考虑到特征值可能为复数，因此将所有特征值均当成复数处理。
此函数中，QR分解部分用子函数void QR_decompositionMk()实现。
这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。
5）计算特征向量得到特征值后，计算实特征值相应的特征向量。
这里判断所得特征值的虚数部分是否为零。
求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。
因此先初始化矩阵Array，计算(A-λI)，再求解方程组。
北航数值分析实验报告
北航数值分析实验报告篇一：北航数值分析报告第一大题《数值分析》计算实习报告第一大题学号：DY1305姓名：指导老师：一、题目要求已知501*501阶的带状矩阵A，其特征值满足?1?2...?501。
试求：1、?1，?501和?s的值；2、A的与数?k??1?k?501??140最接近的特征值?ik（k=1，2，...，39）；3、A的（谱范数）条件数c nd(A)2和行列式de tA。
二、算法设计方案题目所给的矩阵阶数过大，必须经过去零压缩后进行存储和运算，本算法中压缩后的矩阵A1如下所示。
?0?0?A1??a1??b??c0b a2bcc bb c............c bb ccb a500b0a 3...a499c?b??a501??0?0??由矩阵A的特征值满足的条件可知?1与?501之间必有一个最大，则采用幂法求出的一个特征值必为其中的一个：当所求得的特征值为正数，则为?501；否则为?1。
在求得?1与?501其中的一个后，采用带位移的幂法则可求出它们中的另一个，且位移量即为先求出的特征值的值。
用反幂法求得的特征值必为?s。
由条件数的性质可得，c nd(A)2为模最大的特征值与模最小的特征值之比的模，因此，求出?1，?501和?s的值后，则可以求得c nd(A)2。
数值分析课后习题答案
0 1
0 10 1 1 0 0 0 1
0 0 12 1 1 2 0 0 0
1 2
0 0 0 1 1 0
1 2
1 2
1 2
1
0 0 0 1 0
1 2
1 2
0
1 2
1 2
0
0
0
341 1 1
2-5.对矩阵A进行LDLT分解和GGT分解，并求解方程组
Ax=b，其中
16 4 8
1
A 4 5 4 , b 2
8 4 22
3
解
16 A 4
4 5
84
44 11
2-3(1).对矩阵A进行LU分解,并求解方程组Ax=b,其中
2 1 1 A1 3 2
4 ,b6
1 2 2
5
解
2 A 1
1 3
1 2
2 11
22
1
5 2
1
3 21来自，所以 A12
1
2 1 1
5 3
2-2(1).用列主元Gauss消元法解方程组
3 2 6x1 4 10 7 0x2 7 5 1 5x3 6
解
3 2 6 4 10 7 0 7 10 7 0 7
r1r2
消元
10 7 0 7 3 2 6 4 0 0.1 6 6.1
r=0.5101-n/3.162…<0.5101-n/3<0.01% 因此只需n=5.即取101/2=3.1623
北航研究生数值分析试题
∗⎞ ⎟的 A1 ⎠
矩阵。
三、（12 分）试用高斯列主元素法求解线性方程组
⎡ 1 3 −2 −4 ⎤ ⎡ x1 ⎤ ⎡3 ⎤ ⎢ 2 6 −7 −10 ⎥ ⎢ x ⎥ ⎢ −2 ⎥ ⎢ ⎥⎢ 2⎥ = ⎢ ⎥ ⎢ −1 −1 5 9 ⎥ ⎢ x3 ⎥ ⎢14 ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥ ⎢ x4 ⎦ ⎥ ⎣ −6 ⎦ ⎣ −3 −5 0 15 ⎦ ⎣ 四、（12 分）利用矩阵 A 的三角分解 A = LU 求解下列方程组 ⎛ 1 2 1 ⎞ ⎛ x1 ⎞ ⎛ 0 ⎞ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ 2 2 3 ⎟ ⎜ x2 ⎟ = ⎜ 3 ⎟ ⎜ −1 −3 0 ⎟ ⎜ x ⎟ ⎜ 2 ⎟ ⎝ ⎠⎝ 3 ⎠ ⎝ ⎠
第一章
1、近似数 x = 0.231 关于真值 x = 0.229 有（ （1）1；（2）2；（3）3；（4）4。
∗
绪论
一、选择题（四个选项中仅有一项符合题目要求，每小题 3 分，共计 15 分) ）位有效数字。
2、取 3 ≈ 1.732 计算 x = ( 3 − 1) ，下列方法中哪种最好？（
4
）
Ax
∞和
A ∞ 的值分别为（
）
3
(1) 8 , 8 ；
(2) 8 , 7 ；
(3) 8 , 6 ；
(4) 7 , 7 。
5 、若解线性代数方程组的 Gauss 部分选主元方法第二步得到的系数矩阵的第三列向量为
(2
6 3 2 −5 4 2 ) ，则第三步主行是（
T
） (4) 第 6 行。
(1) 第 2 行；
1 − cos x , sin x
x ≠ 0且 x << 1 ；
（2）
1 1− x , − 1+ 2x 1+ x
数值分析第四版课后答案答案第八章
第八章 常微分方程初值问题数值解法1、解：欧拉法公式为221(,)(100),0,1,2+=+=++=n n n n n n n y y hf x y y h x y n代00y =入上式，计算结果为 123(0.1)0.0,(0.2)0.0010,(0.3)0.00501≈=≈=≈=y y y y y y2、解：改进的欧拉法为1112[(,)(,(,))]n n n n n n n n y y h f x y f x y hf x y ++=+++将2(,)=+-f x y x x y 代入上式，得2111111221n n n n n n h hh x x x x y h y +++)+[(-)(+)+(+)]=(-+ 同理，梯形法公式为211122[(1)(1)]-+++++=++++h h n nn n n n h h y y x x x x 将00,0.1y h ==代入上二式，，计算结果见表9—5表 9—5可见梯形方法比改进的欧拉法精确。
3、证明：梯形公式为111[(,)(,)]2n n n n n n hy y f x y f x y +++=++代(,)f x y y =-入上式，得11[]2++=+--n n n n hy y y y解得21110222()()()222n n n n h h h y y y y h h h++----===⋯=+++ 因为01y =，故2()2nn h y h-=+ 对0x∀>，以h 为步长经n 步运算可求得()y x 的近似值n y ，故,,xx nh n h==代入上式有2()2x hn hy h-=+22220000222lim lim()lim(1)lim[(1)]222x x h h xx h h h h hn h h h h h h h y e h h h+-+→→→→-==-=-=+++4、解：令2()xt y x e dt =⎰，则有初值问题2',(0)0x y e y ==对上述问题应用欧拉法，取h=0.5,计算公式为210.5,0,1,2,3n x n n y y e n +=+=由0(0)0,y y ==得1234(0.5)0.5,(1.0) 1.142012708(1.5) 2.501153623,(2.0)7.245021541≈=≈=≈=≈=y y y y y y y y5、解： 四阶经典龙格-库塔方法计算公式见式（9.7）。
北航2010-2015年研究生数值分析报告期末模拟试卷与真题
北航2010-2015年研究生数值分析报告期末模拟试卷与真题数值分析模拟卷A一、填空（共30分，每空3分）1 设-=1511A ，则A 的谱半径=)(a ρ______，A 的条件数)(1A cond =________. 2 设 ,2,1,0,,53)(2==+=k kh x x x f k ，则],,[21++n n n x x x f ＝________, ],,[321+++n n n n x x x x f ，＝________.3 设≤≤-++≤≤+=21,1210,)(2323x cx bx x x x x x S ，是以0，1，2为节点的三次样条函数，则b=________,c=________.4 设∞=0)]([k k x q 是区间［0，1］上权函数为x x =)(ρ的最高项系数为1的正交多项式族，其中1)(0=x q ，则?=10)(dx x xq k ________，=)(2x q ________.5 设=11001a a a a A ，当∈a ________时，必有分解式，其中L 为下三角阵，当其对角线元素)3,2,1(=i L ii 满足条件________时，这种分解是唯一的.二、（14分）设49,1,41,)(21023====x x x x x f , （1）试求)(x f 在]49,41[上的三次Hermite 插值多项式)(x H 使满足2,1,0),()(==i x f x H i i ，)()(11x f x H '='.（2）写出余项)()()(x H x f x R -=的表达式.三、（14分）设有解方程0cos 2312=+-x x 的迭代公式为n n x x cos 3241+=+，（1）证明R x ∈?0均有?∞→=x x n x lim （?x 为方程的根）；（2）取40=x ，用此迭代法求方程根的近似值，误差不超过，列出各次迭代值；（3）此迭代的收敛阶是多少？证明你的结论.四、(16分) 试确定常数A ，B ，C 和，使得数值积分公式有尽可能高的代数精度. 试问所得的数值积分公式代数精度是多少？它是否为Gauss 型的？五、（15分）设有常微分方程的初值问题=='00)(),(y x y y x f y ，试用Taylor 展开原理构造形如)()(11011--++++=n n n n n f f h y y y ββα的方法，使其具有二阶精度，并推导其局部截断误差主项.六、（15分）已知方程组b Ax ＝，其中= ??=21,13.021b A ，（1）试讨论用Jacobi 迭代法和Gauss-Seidel 迭代法求解此方程组的收敛性.（2）若有迭代公式)()()()1(b Ax a x x k k k ++=+，试确定一个的取值围，在这个围任取一个值均能使该迭代公式收敛.七、（8分）方程组，其中，A 是对称的且非奇异.设A 有误差，则原方程组变化为，其中为解的误差向量，试证明 .其中1λ和2λ分别为A 的按模最大和最小的特征值.数值分析模拟卷B填空题（每空2分，共30分）1. 近似数231.0=*x 关于真值229.0=x 有____________位有效数字；2. 设)(x f 可微，求方程)(x f x =根的牛顿迭代格式是_______________________________________________；3. 对1)(3++=x x x f ，差商=]3,2,1,0[f _________________；=]4,3,2,1,0[f ________；4. 已知???? ??-='-=1223,)3,2(A x ，则=∞||||Ax ________________，=)(1A Cond ______________________ ；5. 用二分法求方程01)(3=-+=x x x f 在区间[0,1]的根，进行一步后根所在区间为_________，进行二步后根所在区间为_________________；6. 求解线性方程组=+=+04511532121x x x x 的高斯—赛德尔迭代格式为_______________________________________；该迭代格式迭代矩阵的谱半径=)(G ρ_______________；7. 为使两点数值求积公式：?-+≈111100)()()(x f x f dx x f ωω具有最高的代数精确度，其求积节点应为=0x _____ , =1x _____,==10ωω__________.8. 求积公式)]2()1([23)(30f f dx x f +≈?是否是插值型的__________，其代数精度为___________。
北航数值分析计算实习报告一
北京航空航天大学《数值分析》计算实习报告第一大题学院：自动化科学与电气工程学院专业：控制科学与工程学生姓名：学号：教师：电话：完成日期：2015年11月6日北京航空航天大学Beijing University of Aeronautics and Astronautics实习题目：第一题 设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。
矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1.求1λ，501λ和s λ的值。
2.求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。
3.求A 的(谱范数)条件数2)A (cond 和行列式detA 。
说明：1.在所用的算法中，凡是要给出精度水平ε的，都取12-10=ε。
2.选择算法时，应使矩阵A 的所有零元素都不储存。
3.打印以下内容： （1）全部源程序；（2）特征值），，39,...,2,1(,s 5011=k k i λλλλ以及A det ,)A (cond 2的值。
4.采用e 型输出实型数，并且至少显示12位有效数字。
一、算法设计方案1、求1λ，501λ和s λ的值。
由于||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤，可知绝对值最大特征值必为1λ和501λ其中之一，故可用幂法求出绝对值最大的特征值λ，如果λ=0，则1λ=λ，否则501λ=λ。
将矩阵A 进行一下平移：I -A A'λ= （1）对'A 用幂法求出其绝对值最大的特征值'λ，则A 的另一端点特征值1λ或501λ为'λ+λ。
北航数值分析计算实习题目一 幂法 反幂法 求矩阵特征值
《数值分析》计算实习题目第一题：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）作为初始向量进行迭代，可得出以上结果。
北航研究生数值分析编程大作业1
数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。
这样所需要的存储单元数大大减少，从而极大提高了运算效率。
2、利用幂法求出5011λλ，幂法迭代格式：011111111nT k k k k k k kk T k k k u R u u y u u Ay y u ηηβ--------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。
首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。
3、利用反幂法求出iks λλ,反幂法迭代格式： 011111111nTk k k k k k kk T k k k u R u u y u Au y y u ηηβ--------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。
每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。
北航数值分析B自测题1-3章
D.2
10.若f (x)在区间[a, b]内单调有限，二分k 次后区间记为[ak , bk ],且每次 取xk+1 = ak 近似代替精确解x∗ ,则最小的绝对误差限为（ ） A.|xk+1 − x∗ | ≤ C.|xk+1 − x∗ | ≤ 二、 填空题 1.设π 的近似数π ∗ 有4位有效数字，则其相对误差限为 √ 2. x∗ 的相对误差约是x∗ 的相对误差的 倍。
(2)a=6,b=2; (3)a=2,b=3; (4)a=-1,b=2.
3.设矩阵A ∈ Rn∗n , Q ∈ Rn∗n ,且QT Q=E，则下列关系式不成立的是（ ） (1) A
2
= AQ ; (2) QA
2 F
= A
F
; (3) Qx
2
= x , 其中x ∈ Rn ;
2
(4)cond∞ ( ) = cond∞ ( Q).
8
3 −1 4 1 4.设矩阵A=−1 2 −2 ，x=−1 ，则 Ax 2 −3 −2 1 （ ） (1)8,8; (2)8,7; (3)8,6; (4)7,7.
A. 7 3
B.− 7 3
7 C. 6
D.− 7 6
8.追赶法适用于求解（ ）线性方程组. A.上三角 B.下三角 C.对角 D.三对角 二.填空 题：
2 1 0 1.设A=1 2 a，为使A 可分解为A=LLT ，其中L是 对角元素为正的下 0 a 2 三角矩阵 范围是 ，则a的取值 .
1 (s∗ ) s∗
≈
27 110×80
= 0.31% = max 1 ≤ j ≤ n{1, 2, 4} = 4;
北京航空航天大学《工科数学分析》考试试题及参考答案(2012-2013第一学期)
f x e
'
e
cos x ln sin x
cos 2 x sin x cos x sin x ln sin x . sin x
dy dy dx cos t t sin t 4）解： . dx cos t t sin t dt dt
m 满足什么条件，函数在 x 0 可导.
2. 证明下面问题（10 分） 设 s 0, x1 0, xn1
1 s x , 证明数列 xn 单调有界，且极限为 s . n 2 x n
1 , 用 Cauchy 收敛定理证明 xn 收敛. 2n
5.
1) 用反证法证明. 假设存在 q a, b , g q 0 . 则根据拉格朗日中值定理
' g a g q g ' x1 a q 0 得到 g x1 0, x1 a, q
g b g q g ' x2 b q 0 得到 g ' x2 0, x2 q , b
7.
(10 分）证明下面问题 设 f x 定义在 a, b 上. 如果对 a, b 内任何收敛的点列 xn 都有 lim f xn 存在, 则
n
f 在 a, b 上一致连续.
8. (10 分）附加题 (下面两个题目任选其一) 1) 设函数 f
n 1 2 n cos x Cn cos 2 x 1 Cn cos n x , x Cn n1
二、第一次考试题目及答案
1. 计算下面各题（满分 40 分，每个题目 5 分） 1) 2) 计算极限 lim
x 0
北航数值分析报告第一次大作业(幂法反幂法)
一、问题分析与算法描述1. 问题的提出：〔1〕用幂法、反幂法求矩阵的按摸最大和最小特征值，并求出相应的特征向量。
其中要求：迭代精度达到。
〔2〕用带双步位移的QR法求上述的全部特征值，并求出每一个实特征值相应的特征向量。
2. 算法的描述：(1) 幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。
其迭代格式为：终止迭代的控制选用。
幂法的使用条件为实矩阵A具有n个线性无关的特征向量，其相应的特征值满足不等式或幂法收敛速度与比值或有关，比值越小，收敛速度越快。
(2) 反幂法反幂法用于计算实矩阵A按摸最小的特征值，其迭代格式为：每迭代一次都要求解一次线性方程组。
当k足够大时，，可近似的作为矩阵A的属于的特征向量。
比值越小，收敛的越快。
反幂法要求矩阵A非奇异。
(3) 带双步位移的QR分解法QR方法适用于计算一般实矩阵的全部特征值，尤其适用于计算中小型实矩阵的全部特征值。
本算例中采用带双步位移的QR方法，可加速收敛，其迭代格式为：二、计算结果与分析1. 计算结果：(1) 幂法：初始条件：最大迭代次数L=1000；向量计算结果：第1次迭代结果：最大特征值：0.00000e+000第2次迭代结果：最大特征值：2.48910e+000 相对误差：1.00000e+000 第3次迭代结果：最大特征值：1.67719e+000 相对误差：第4次迭代结果：最大特征值：-2.10960e+000 相对误差：1.79503e+000 第5次迭代结果：最大特征值：-6.13203e-001 相对误差：2.44030e+000 ……第794次迭代结果：最大特征值：-1.97638e+000 相对误差：最大特征值：-1.97638e+000 相对误差：********************最终迭代结果***************特征值：-1.97638e+000 相对误差：迭代次数：795(2) 反幂法：初始条件：最大迭代次数L=1000；向量运行结果：第1次迭代结果：最大特征值：1.07542e+000第2次迭代结果：最大特征值：-3.66550e+000 相对误差：1.29339e+000 第3次迭代结果：最大特征值：1.22709e+001 相对误差：1.29871e+000 第4次迭代结果：最大特征值：-1.03421e+000 相对误差：1.28650e+001 第5次迭代结果：最大特征值：相对误差：……第995次迭代结果：最大特征值：相对误差：第996次迭代结果：最大特征值：相对误差：最大特征值：相对误差：第998次迭代结果：最大特征值：相对误差：第999次迭代结果：最大特征值：相对误差：第1000次迭代结果：最大特征值：相对误差：******************************超过最大设定迭代次数，迭代失败！(3) 带双步位移的QR法：初始条件：最大迭代次数L=1000；向量运行结果：全部特征值：特征向量〔经谱X数归一化〕：实特征值对应特征向量：-0.062705 -0.022368 0.304372 0.064466 0.521833 -0.157024 0.136942 -0.218108 0.250264 -0.043064 -0.228688 -0.184632 -0.072871 0.124721 0.029070 0.102566 -0.136358 0.167727 0.085747 0.546165 实特征值对应特征向量：-0.018001 0.019652 0.273447 0.070528 0.274896 -0.144015 0.048385 0.376439 -0.583051 -0.054008 -0.168682 -0.113430 -0.034709 0.009204 0.472291 0.125664 -0.190617 0.113145 0.046278 0.059871 实特征值对应特征向量：0.106861 0.087709 -0.024967 -0.020897 0.064302 0.034047 0.535143 0.046383 0.028832 0.003479-0.097276 -0.383801 0.089445 -0.039560 -0.036928 -0.021330 0.014811 0.705836 -0.108904 0.082022 实特征值对应特征向量：-0.055201 0.003399 0.242191 0.102847 0.372470 -0.372826 0.113953 0.240659 -0.310401 -0.076590 -0.244632 -0.192549 -0.077259 0.263328 0.201662 0.154166 -0.407814 0.186782 0.094649 0.173302 实特征值对应特征向量：0.427828 -0.546801 0.007822 -0.382580 0.025199 0.012788 0.033241 0.005389 -0.004065 0.043524 -0.032112 -0.044233 0.135395 -0.006564 0.001214 0.020165 0.011678 0.050001 -0.585765 0.013115 实特征值对应特征向量：0.236032 -0.139250 -0.008143 0.638527 -0.009049 -0.002911 -0.001307 0.003054 0.006515 -0.030134 0.012712 0.011368 -0.018792 -0.001753 -0.005749 -0.014290 -0.005292 -0.014591 0.717590 0.001369 实特征值对应特征向量：-0.227404 -0.048154 0.022615 0.297305 0.070372 0.039927 0.078503 0.015822 -0.012182 0.605334 -0.083616 -0.106270 -0.573963 -0.019907 0.003839 0.051362 0.036567 0.115613 0.332707 0.036954 实特征值对应特征向量：-0.027768 -0.051081 -0.159642 -0.054573 -0.084441 0.118378 0.029553 0.211088 0.203867 0.0486272. 结果分析以上三种方法中，幂法计算共进展了795次迭代才达到收敛，计算量较大，收敛性不好；反幂法计算结果未能收敛，通过进一步分析发现，这是因为反幂法迭代程序未考虑按模最小特征值为复数的情况，造成迭代失败。
北航数值分析A大作业
一、算法设计方案1、解非线性方程组将各拟合节点（x i ，y j ）分别带入非线性方程组，求出与(,)i i x y 相对应的数组te[i][j]，ue[i][j]，求解非线性方程组选择Newton 迭代法，迭代过程中需要求解线性方程组，选择选主元的Doolittle 分解法。
2、二元二次分偏插值对数表z(t,u)进行分片二次代数插值，求得对应（t ij ，u ij ）处的值，即为),(j i y x f 的值。
根据给定的数表，可将整个插值区域分成 16 个小 的区域，故先判断t ij , u ij 所在，的区域，再作此区域的插值，计算 z ij ，相应的Lagrange 形式的插值多项式为：112211(,)()()(,)m n krkrk m r n p t u l t l u f t u ++=-=-=∑∑其中11()m wk w m k ww kt t l t t t +=-≠-=-∏ (k=m-1, m, m+1) 11()n wr w n r ww ry y l u y y +=-≠-=-∏ (r=n-1, n, n+1)3、曲面拟合从k=1开始逐渐增大k 的值，使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，当710-<σ时结束计算。
拟合基函数φr (x)ψs (y)选择为φr (x)=x r ，ψs (y)=y s 。
拟合系数矩阵c 通过连续两次解线性方程组求得。
[]rsc *=C ，11()()T T T --=C B B B UG G G其中0011101011[()]1kk r i k x x x x x x x ϕ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦B ，0011101011[()]1kk s j k y y y y G y y y ψ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦[(,)]i j f x y =U4、观察比较计算)5,,2,1,8,,2,1)(,(),,(****⋅⋅⋅=⋅⋅⋅=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。
北航数值分析-实习作业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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
北京航空航天大学数值分析大作业八学院名称自动化专业方向控制工程学号学生姓名许阳教师孙玉泉日期2014 年11月26 日一．题目关于x , y , t , u , v , w 的方程组(A.3)⎪⎪⎩⎪⎪⎨⎧=-+++=-+++=-+++=-+++79.0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。
表A-1 二维数表1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式∑∑===k i kj s r rs y x c y x p 00),(要求p (x , y )以最小的k 值达到以下的精度∑∑==-≤-=100207210)],(),([i j i i i i y x p y x f σ其中j y i x i i 05.05.0,08.0+==。
2. 计算),(),,(****j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。
二．算法设计（一）总体思路1.题目要求∑∑===ki kj s r rs y x c y x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。
),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。
2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。
（二）算法实现1. ),(i i y x 与),(i i y x f 的数表的获得对区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的f (x , y )值可由方程组及二维数表得到。
将区域D 上的),(i i y x 分别回代于方程组(A.3)，成为关于t,u,v,w 的4元非线性方程组，解出每个),(j i y x 对应的t,u 。
再通过表A-1进行插值近似，得到相应的z 值。
对应的z 即为D 区域上),(j i y x 对应的),(j i y x f 。
从而得到),(j i y x 与),(i i y x f 的数表。
(1) 4元非线性方程组求解),(i i y x 代入(A.3)后，原方程组变为关于t,u,v,w 的4元非线性方程组。
观察到方程组中方程形式较为简单，易于对变量t,u,v,w 求偏导数，故而选用Newton 法对方程组求解。
计算方程组矩阵为：⎪⎪⎪⎪⎪⎭⎫⎝⎛--+++--+++--+++--+++=79.0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0),,,(i i i i y w v u t x w v u t y w v u t x w v u t w v u t F计算方程组偏导数矩阵为：⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--=w v u tw v u t F cos 15.011sin 15.011cos 5.01111sin 5.0),,,('迭代公式为：⎪⎪⎪⎪⎪⎭⎫⎝⎛∆∆∆∆+⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛++++)()()()()()()()()1()1()1()1(k k k k k k k k k k k k w v u t w v u t w v u t ，k=0,1,2,…,n 其中⎪⎪⎪⎪⎪⎭⎫ ⎝⎛∆∆∆∆)()()()(k k k k w v u t 为线性方程组),,,(),,,(')()()()()()()()()()()()(k k k k k k k k k k k k w v u t F w v u t w v u t F -=⎪⎪⎪⎪⎪⎭⎫⎝⎛∆∆∆∆的解。
取ε≤∆∆∆∆|)||,||,||,m ax (|/|)||,||,||,m ax (|)()()()()()()()(k k k k k k k k w v u t w v u t 为迭代终止条件。
由表A-1观察到t,u 基本在（0,2）上，于是选取)1,1,1,1(),,,)0()0()0()0(=w v u t （为迭代初值。
通过以上方法求得与),(j i y x 对应的),t (ij ij u 。
(2) 分片二元双二次代数插值为保证代数插值的收敛性，应采用分片低次插值。
故此使用分片双二次代数插值。
)5,...,1,0;5,...,1,0(4.0,2.0====n m n u m t n m给定),t (ij ij u 如满足如下关系式：1.01.0+≤<-m ij m t t t ，41<<m2.02.0+≤<-n ij n u u u ，41<<n则选择)2,1,;2,1,)(,(++=++=n n n r m m m k u t r k 为插值节点，相应插值多项式为∑∑+=+==2222),()()(),(m m k n n r r k ij r ij k ij ij u t z u l t l u t p其中∏+≠=--=2)(m k q m q qkq ij ij k t t t t t l )21,++=m m m k ，（∏+≠=--=2)(n kq n q qk q ij ij r u u u u u l )2,1,++=n n n r （如果7.03.0>≤ij ij t t 或，则上式取m=1或m=4；如果6.0≤ij u 或4.1>ij u ，则取n=1或n=4。
得到),(22ij ij u t p 表达式后，直接带入),t (ij ij u ，得到的值即为与),(j i y x 对应的),(j i y x f 。
2. 乘积型最小二乘曲面拟合使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。
UG B G G C B B T T T =)()(3.计算),(),,(****j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值),(**j i y x f 的计算与),(j i y x f 相同。
将),(**j i y x 代入原方程组，求解响应),(**ij ij u t 进行分片双二次插值求得),(**j i y x f 。
),(**j i y x p 的计算则可以直接将),(**j i y x 代入所求p(x,y)。
三．Matlab源程序及结果牛顿法解非线性方程组子程序：function [t,u]=newt(x,y)t=1; u=1; v=1; w=1;ep=1e-12;k=1;N=100;while(k<N)F(1,1)=0.5*cos(t)+u+v+w-x-2.67;F(2,1)=t+0.5*sin(u)+v+w-y-1.07;F(3,1)=0.5*t+u+cos(v)+w-x-3.74;F(4,1)=t+0.5*u+v+sin(w)-y-0.79;dF=[-0.5*sin(t) 1 1 1;1 0.5*cos(u) 1 1;0.5 1 -sin(v) 1;1 0.5 1 cos(w)]; deltaX=ones(4,1);deltaX=dF^-1*(-1)*F;if max(abs(deltaX))/abs(x)<ept=t+deltaX(1,1);u=u+deltaX(2,1);v=v+deltaX(3,1);w=w+deltaX(4,1);break;endt=t+deltaX(1,1);u=u+deltaX(2,1);v=v+deltaX(3,1);w=w+deltaX(4,1);k=k+1;end分片双二次代数插值子程序：function f=p22(t,u)z=[-0.5 -0.34 0.14 0.94 2.06 3.5;-0.42 -0.5 -0.26 0.3 1.18 2.38;-0.18 -0.5 -0.5 -0.18 0.46 1.42;0.22 -0.34 -0.58 -0.5 -0.1 0.62;0.78 -0.02 -0.5 -0.66 -0.5 -0.02;1.5 0.46 -0.26 -0.66 -0.74 -0.5];if (t<=0.3)i=1;endif (t>0.3&&t<=0.5)i=2;endif (t>0.5&&t<=0.7)i=3;endif (t>0.7)i=4;endif u<=0.6j=1;endif u>0.6&&u<=1j=2;endif u>1&&u<=1.4j=3;endif u>1.4j=4;endfor k=i:i+2num=1;for a=i:i+2if a~=knum=num*(t-0.2*(a-1))/(0.2*(k-1)-0.2*(a-1));endendl(k)=num;endfor r=j:j+2num=1;for a=j:j+2if a~=rnum=num*(u-0.4*(a-1))/(0.4*(r-1)-0.4*(a-1));endendm(r)=num;endsum=0;for k=i:i+2for r=j:j+2sum=sum+l(k)*m(r)*z(k,r);endendf=sum;最小二乘曲面拟合子程序：function [C,k,sigma]=fitxy(A)k=1;N=10;while k<NB=ones(11,k+1);G=ones(21,k+1);for i=1:kfor l=1:11B(l,i+1)=(0.08*(l-1))^i;endfor m=1:21G(m,i+1)=(0.5+0.05*(m-1))^i;endendU=zeros(11,21);for i=1:11for j=1:21U(i,j)=A((i-1)*21+j,3);endendC=(B'*B)^-1*B'*U*G*(G'*G)^-1;sigma=0;for i=1:11for j=1:21p=0;for r=0:kfor s=0:kp=p+C(r+1,s+1)*(0.08*(i-1))^r*(0.5+(0.05*(j-1)))^s;endendsigma=sigma+(U(i,j)-p)^2;endEndksigmaif sigma<=1e-7break;endk=k+1;end主程序：clc;clear;A1=zeros(11*21,3);for i=1:11x(i)=0.08*(i-1);for j=1:21y(j)=0.5+0.05*(j-1);[t(i,j),u(i,j)]=newt(x(i),y(j));A1((i-1)*21+j,1)=x(i);A1((i-1)*21+j,2)=y(j);A1((i-1)*21+j,3)=p22(t(i,j),u(i,j));endend[C,k,sigma]=fitxy(A1);A2=zeros(40,4);for i=1:8x1(i)=0.1*i;for j=1:5y1(j)=0.5+0.2*j;[t1(i,j),u1(i,j)]=newt(x1(i),y1(j));A2((i-1)*5+j,1)=x1(i);A2((i-1)*5+j,2)=y1(j);A2((i-1)*5+j,3)=p22(t1(i,j),u1(i,j));endendfor i=1:8for j=1:5p=0;for r=0:kfor s=0:kp=p+C(r+1,s+1)*(0.1*i)^r*(0.5+(0.2*j))^s;endendA2((i-1)*5+j,4)=p;endendA1A2程序结果：数表：)),(,,i i i i y x f y x （k和 :数表)),(),,(,,(******j i j i j i y x p y x f y x ：。