北航数值分析大作业第二题(f o r t r a n)
-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN
“数值分析“计算实习大作业第二题
——SY1415215
孔维鹏一、计算说明
本程序采用带双步位移的QR方法求解矩阵A的所有特征值，然后采用反幂法求解矩阵A的实特征值对应的特征向量。
在采用带双步位移的QR方法求解特征值时，对教材上所提供的具体算法作稍微的改动，以简化程序，具体算法如下所示：
1、计算出A拟上三角化后的矩阵，给定精度水平和最大迭代次数L；
2、记，令k=1，m=n；
3、如果，则可直接计算出最后1或2个特征值，转8，否则转4；
4、如果，则可得一个特征值，置m=m-1；转3，否则转5；
5、如果，则可得两个特征值，置m=m-2；转3，否则转6；
6、记，计算
7、k=k+1,转3
8、A的全部特征值已经求出，停止计算。
二、计算源程序（FORTRAN）
PROGRAM SY1415215_2
PARAMETER (N=10)
DIMENSION A(N,N),A1(N,N),A2(N,N),C(2,N),Q(N,N),R(N,N),CR(N),CM(N)!C为存储特征值的数组,1为实部，为虚部
REAL(8) A,A1,A2,C,Q,R,CM
E=1E-12 !精度水平
L=1000 !迭代最大次数
OPEN(1,FILE='数值分析大作业第二题计算结果.TXT')
DO I=1,N
DO J=1,N
IF(I==J) THEN
A(I,J)=*COS(I+*J)
ELSE
A(I,J)=SIN*I+*J)
ENDIF
ENDDO
ENDDO
A1=A
WRITE(*,"('矩阵A为：')")
WRITE(1,"('矩阵A为：')")
DO I=1,N
DO J=1,N
WRITE(*,"(2X,,2X,\)") A(I,J)
WRITE(1,"(2X,,2X,\)") A(I,J)
ENDDO
WRITE(*,"(' ')")
WRITE(1,"(' ')")
ENDDO
!使用矩阵的拟上三角化的算法将矩阵A化为拟上三角矩阵A（n-1）CALL HESSENBERG(A,N)
WRITE(*,"('拟上三角化后矩阵A（n-1）为：')")
WRITE(1,"('拟上三角化后矩阵A（n-1）为：')")
DO I=1,N
DO J=1,N
WRITE(*,"(2X,,2X,\)") A(I,J)
WRITE(1,"(2X,,2X,\)") A(I,J)
ENDDO
WRITE(*,"('')")
WRITE(1,"('')")
ENDDO
!计算对矩阵A（n-1）实行QR方法迭代结束后所得矩阵
A2=A
CALL QRD(A2,N,Q,R)
WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得Q为：')") DO I=1,N
DO J=1,N
WRITE(*,"(2X,,2X,\)") Q(I,J)
WRITE(1,"(2X,,2X,\)") Q(I,J)
ENDDO
WRITE(*,"('')")
WRITE(1,"('')")
ENDDO
WRITE(*,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") WRITE(1,"('对矩阵A（n-1）实行QR方法迭代结束后所得R为：')") DO I=1,N
DO J=1,N
WRITE(*,"(2X,,2X,\)") R(I,J)
WRITE(1,"(2X,,2X,\)") R(I,J)
ENDDO
WRITE(*,"('')")
WRITE(1,"('')")
ENDDO
!使用带双步位移的QR方法求解矩阵A（n-1）的特征值
K=1
M=N
DO WHILE(K<=L)
IF(M<=2) THEN
IF(M==1) THEN
C(1,M)=A(M,M)
ELSE IF(M==2) THEN
CALL CALCUS(A,N,M,C)
ENDIF
EXIT
ELSE IF(ABS(A(M,M-1)) C(1,M)=A(M,M) M=M-1 ELSE IF(ABS(A(M-1,M-2)) CALL CALCUS(A,N,M,C) M=M-2 ELSE CALL CALM(A,M,N) ENDIF K=K+1 ENDDO WRITE(*,"('矩阵A的全部特征值为：')") WRITE(1,"('矩阵A的全部特征值为：')") DO J=1,N WRITE(*,",'+',,'i')") C(1,J),C(2,J) WRITE(1,",'+',,'i')") C(1,J),C(2,J) ENDDO !使用反幂法求解A的相应于实特征值的特征向量J=1 DO I=1,N IF(C(2,I)==0)THEN CR(J)=C(1,I) J=J+1 ENDIF ENDDO JC=J-1 WRITE(*,"('矩阵A的实特征值为：')") 《数值分析B》大作业一 SY1103120 朱舜杰 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（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和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和 最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。 二．源程序 #include 数值分析作业答案 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange插值基底。 （3）用Newton基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为: ， 所以： 所以f(x)的二次插值多项式为： （2）用Lagrange插值基底 Lagrange插值多项式为： 所以f(x)的二次插值多项式为： (3) 用Newton基底: 均差表如下： xk f(xk) 一阶均差二阶均差 1 0 -1 -3 3/2 2 4 7/ 3 5/6 Newton插值多项式为： 所以f(x)的二次插值多项式为： 由以上计算可知，三种方法得到的多项式是相同的。 6、在上给出的等距节点函数表，若用二次插值求ex的近似值，要使截断误差不超过10-6，问使用函数表的步长h应取多少？ 解：以xi-1,xi,xi+1为插值节点多项式的截断误差，则有 式中 令得 插值点个数 是奇数，故实际可采用的函数值表步长 8、，求及。 解：由均差的性质可知，均差与导数有如下关系： 所以有： 15、证明两点三次Hermite插值余项是 并由此求出分段三次Hermite插值的误差限。 证明：利用[xk，xk+1]上两点三次Hermite插值条件 知有二重零点xk和k+1。设 确定函数k(x): 当或xk+1时k(x)取任何有限值均可； 当时，，构造关于变量t的函数 显然有 在[xk,x][x,xk+1]上对g(x)使用Rolle定理，存在及使得 在，，上对使用Rolle定理，存在，和使得 再依次对和使用Rolle定理，知至少存在使得 而，将代入，得到 推导过程表明依赖于及x 综合以上过程有： 确定误差限： 记为f(x)在[a,b]上基于等距节点的分段三次Hermite插值函数。在区间[xk，xk+1]上有 而最值 进而得误差估计： 16、求一个次数不高于4次的多项式，使它满足，，。 《数值分析》计算实习题目 第一题： 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 目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。 这样即求出所有A所有实特征值对应的一个特征向量。 #include 数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： x y当作已知量代入题目给定的非线性方程组，求（1）解非线性方程组。将给定的(,) i i 得与(,)i i x y 相对应的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解 (,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 （二）具体算法设计： （1）解非线性方程组 牛顿法解方程组()0F x =的解* x ，可采用如下算法： 1）在* x 附近选取(0) x D ∈，给定精度水平0ε>和最大迭代次数M 。 2）对于0,1, k M =执行 ① 计算() ()k F x 和()()k F x '。 ② 求解关于() k x ?的线性方程组 () ()()()()k k k F x x F x '?=- ③ 若() () k k x x ε∞∞ ?≤，则取*()k x x ≈，并停止计算；否则转④。 ④ 计算(1) ()()k k k x x x +=+?。 ⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。 （2）分片双二次插值 给定已知数表以及需要插值的节点，进行分片二次插值的算法： 设已知数表中的点为： 00(0,1,,) (0,1,,)i j x x ih i n y y j j m τ=+=???=+=?? ，需要插值的节点为(,)x y 。 1) 根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+ 或12 n h x x ->-，插值节点对应取1i =或1i n =-， 数值分析大作业（2013年5月） 金洋洋（12721512），机自系 1.下列各数都是经过四舍五入得到的近似值，试分别指出它 们的绝对误差限, 相对误差限和有效数字的位数。 X1 =5.420, x 2 =0.5420, x 3=0.00542, x 4 =6000, x 5=50.610? 解：根据定义：如果*x 的绝对误差限 不超过x 的某个数位的半个单位，则从*x 的首位非零数字到该位都是有效数字。 显然根据四舍五入原则得到的近视值，全部都是有效数字。 因而在这里有：n1=4, n2=4, n3=3, n4=4, n5=1 （n 表示x 有效数字的位数） 对x1：有a1=5, m1=1 （其中a1表示x 的首位非零数字，m1表示x1的整数位数） 所以有绝对误差限 143 11 (1)101022 x ε--≤ ?=? 相对误差限 31() 0.510(1)0.00923%5.4201 r x x x εε-?= == 对x2：有a2=5, m2=0 所以有绝对误差限 044 11 (2)101022 x ε--≤ ?=? 相对误差限 42() 0.510(2)0.00923%0.54202 r x x x εε-?= == 对x3：有a3=5, m3=-2 所以有绝对误差限 235 11 (3)101022 x ε---≤ ?=? 相对误差限 53() 0.510(3)0.0923%0.005423 r x x x εε-?= == 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022 x ε-≤?= 相对误差限 4() 0.5 (4)0.0083%6000 4 r x x x εε= = = 对x5：有a5=6, m5=5 所以有绝对误差限 514 11(5)101022 x ε-≤ ?=? 相对误差限 45() 0.510(5)8.3%600005 r x x x εε?= == 《数值分析》计算实习报告 第二题 院系：机械工程及自动化学院 学号： 姓名： 2017年11月 一、题目要求 试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知 a ij ={ sin (0.5i +0.2j ) i ≠j 1.52cos (i +1.2j ) i =j (i,j =1,2, (10) 说明： 1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10?12。 2.打印以下内容： （1）全部源程序； （2）矩阵A 经过拟上三角化后所得的矩阵A (n?1)； （3）对矩阵A (n?1)实行QR 方法迭代结束后所得的矩阵； （4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,?,10)，其中R i =Re(λi )，I i = Im(λi ) 。若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。 3.采用e 型数输出实型数，并且至少显示12位有效数字。 二、算法设计思路和方案 1. 将矩阵A 拟上三角化得到矩阵A (n?1) 为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n?1)，然后用QR 方法计算A (n?1)的全部特征值，而A (n?1)的特征值就是A 的特征值。具体算法如下： 记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2, ,;,1,,)ij a i n j r r n ==+。 对于1,2,,2r n =-执行 （1）若() (2,3,,)r ir a i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。 数值分析第二次大作业 史立峰 SY1505327 一、 方案 （1）利用循环结构将sin(0.50.2)() 1.5cos( 1.2)() {i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的 矩阵A ； （2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。 对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij ,,1,;,,2,1) ( +==。 对于2,,2,1-=n r 执行 1. 若 ()n r r i a r ir ,,3,2) ( ++=全为零，则令A(r+1) =A(r),转5；否则转2。 2. 计算 () ∑+== n r i r ir r a d 1 2 )( ()( )r r r r r r r r r r d c a d a c ==-=++则取,0sgn ) (,1)(,1若 )(,12r r r r r r a c c h +-= 3. 令 () n T r nr r r r r r r r r R a a c a u ∈-=++) ()(,2)(,1,,,,0,,0 。 4. 计算 r r T r r h u A p /)(= r r r r h u A q /)(= r r T r r h u p t /= r r r r u t q -=ω T r r T r r r r p u u A A --=+ω)()1( 5. 继续。 （3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下： 1. 给定精度水平0>ε和迭代最大次数L 。 2. 记n n ij n a A A ?-==][) 1()1()1(，令n m k ==,1。 题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵 20010-1-24A=0-2131 43 1?? ? ? ? ??? 为例进行具体的求解计算。 一、 算法分析： 一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。 二、 算法设计： 1、 原始实矩阵A R n n ?∈拟上三角化 为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’ 记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于r=1,2,…,n -2执行 (1) 若() (2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。 (2) 计算 令 ()2 ()() 1,1,1,sgn(0,sgn()=1) r r r r r r r r r r r c a a a ρ+++=-=，（若则取 (3) 令-0=r n r u u ?? ? ?? ，()()()-1,2,1(,,...,)r r r T n r r r r r r nr u a c a a ρ++=- (4) 计算 (r)(r)(r)T n-r r+1,r r r+2,r nr r T n-r T n-r n-r n-r n-r r+1r 1u = (a -c ,a ,...,a )ρ I H =I -2uu =H H =I -2u u A =HA H ?? ? ?? (5) 继续 算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。 2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值 记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵， 《数值分析A》 一、算法设计方案 整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。 1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。 2.对经过拟上三角化处理的矩阵进行QR分解。 3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。 二、源程序代码 #include 北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期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 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([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. 给定初值0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0) 《数值分析A》计算实习题目第一题 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（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和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有 对角线上元素的乘积。 二．源程序(VS2010环境下，C++语言) #include 数值分析上机作业 第 1 章 1.1计算积分，n=9。（要求计算结果具有6位有效数字） 程序： n=1:19; I=zeros(1,19); I(19)=1/2*((exp(-1)/20)+(1/20)); I(18)=1/2*((exp(-1)/19)+(1/19)); for i=2:10 I(19-i)=1/(20-i)*(1-I(20-i)); end format long disp(I(1:19)) 结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。取6位有效数字可得. 1.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf 命令画出二元函数 z= 的三维图形。 程序： >> x = -10:0.1:10; y = -10:0.1:10; [X,Y] = meshgrid(x,y); Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1); subplot(2,2,1); mesh(X,Y,Z); title('步长0.1') >> x = -10:0.2:10; y = -10:0.2:10; [X,Y] = meshgrid(x,y); Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1); subplot(2,2,1); mesh(X,Y,Z); title('步长 0.2') >>x = -10:0.05:10; y = -10:0.05:10; [X,Y] = meshgrid(x,y); Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1); subplot(2,2,1); mesh(X,Y,Z); title('步长0.05') 第2章 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange 插值基底。 （3）用Newton 基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为:2 210)(x a x a a x P ++=， 所以：64 211111 1111122 2 211 200 -=-==x x x x x x A 3 76144 211111114241 13110111)() ()(22 221120 022 2 22 11 120 00-=-= ---==x x x x x x x x x f x x x f x x x f a 2 3694211111114411 31101111)(1)(1 )(122 221120 02 2 22112 001=--= --==x x x x x x x x f x x f x x f a 6 5654 2 1 1111114 2 1 3 11011111) (1)(1)(122 2 21120 022 11 00 2=--= ---==x x x x x x x f x x f x x f x a 所以f(x)的二次插值多项式为：26 52337)(x x x P ++-= （2）用Lagrange 插值基底 )21)(11() 2)(1())(())(()(2010210-+-+=----=x x x x x x x x x x x l )21)(11() 2)(1())(())(()(2101201------=----=x x x x x x x x x x x l ) 12)(12() 1)(1())(())(()(1202102+-+-=----= x x x x x x x x x x x l 一、算法设计方案 1.使用牛顿迭代法，对原题中给出的i x i 08.0=，j y j 05.05.0+=， (010 ,020i j ≤≤≤≤)的11*21组j i y x ,分别求出原题中方程组的一组解，于是得到一组和i i y x ,对应的j i t u ,。 2.对于已求出的j i t u ,，使用分片二次代数插值法对原题中关于u t z ,,的数表进行插值得到 ij z 。于是产生了z=f(x,y)的11*21个数值解。 3.从k=1开始逐渐增大k 的值，并使用最小二乘法曲面拟合法对z=f(x,y)进行拟合，得到每次的σ,k 。当7 10-<σ时结束计算，输出拟合结果。 4.计算)5,,2,1,8,,2,1)(,(),,(* ***???=???=j i y x p y x f j i j i 的值并输出结果，以观察),(y x p 逼近),(y x f 的效果。其中j y i x j i 2.05.0,1.0* *+==。 二、算法实现方案 1、求(,)f x y ： （1）Newton 法解非线性方程组 0.5cos 2.670.5sin 1.07(1)0.5cos 3.740.5sin 0.79 t u v w x t u v w y t u v w x t u v w y +++-=??+++-=? ? +++-=??+++-=?， 其中，t, u, v ,w 为待求的未知量，x, y 为代入的已知量。 设(,,,)T t u v w ξ=，给定精度水平12110ε-=和最大迭代次数M ，则解该线性方程组的迭代格式为： *(0)(0)(0)(0)(0)(k+1) ()()1()(,,,)()()0,1,T k k k t u v w F F k ξξξ ξξξ-?=?'=-??= ? 在附近选取初值， 迭代终止条件为()(1) () 1/k k k ξξ ξε-∞ ∞ -≤，若k M >时仍未达到迭代精度，则迭代计算失 败。 其中，雅可比矩阵 0.5*cos(t) + u + v + w - x - 2.67t + 0.5*sin(u) + v + w - y - 1.07()0.5*t + u + cos(v) + w - x - 3.74t + 0.5*u + v + sin(w) - y - 0.79F ξ???? ? ?=?????? ， 《数值分析》计算实习作业 （第二题） 算法设计方案： 1、对矩阵A 赋值，取计算精度ε=1×10-12； 2、对矩阵A 进行拟上三角化，得到A (n-1)，并输出A (n-1)； 对矩阵A 的拟上三角化，通过直接调用子函数inftrianglize(A)来实现；拟上三角化得到的矩阵A (n-1)输出至文件solution.txt 中。 3、对A (n-1)进行QR 分解并输出Q 、R 及RQ 矩阵； QR 分解通过直接调用子函数QRdescom(A,Q,R, n)实现。 4、运用QR 方法求所有的特征值，并输出； (1)初始时令m=n ，在m>2的条件下执行； (2)判断如果|A mm-1|<ε，则得到一个特征值，m=m-1，转(4)；否则转(3)； (3)判断如果|A m-1m-2|<ε，则得到两个特征值，m=m-2，转(4)； (4)判断如果m ≤2，转(6)；否则转(5)； (5)执行相似迭代，转(2)； k k T k k k k k k k k k k Q A Q A R Q M I D A D tr A M ==+-=+1)2)det(( (6)求出最后的一个或两个特征值； (7)输出全部的特征值至文件solution.txt 中。 5、输出QR 分解法迭代结束之后的A (n-1)至文件solution.txt 中； 6、通过反幂法求出所有实特征值的特征向量并输出。 首先令B=(A-λi I)，其中λi 是实特征值；反幂法通过调用子函数Bpowmethod(B,x1)实现，最终λi 对应的特征向量就是x1；最后将所有的实特征值的特征向量输出。 北京航空航天大学2020届研究生 《数值分析》实验作业 第九题 院系：xx学院 学号： 姓名： 2020年11月 Q9:方程组A.4 一、 算法设计方案 （一）总体思路 1.题目要求∑∑=== k i k j 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 的数表的获得 对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。再将 ),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。 2.乘积型最小二乘曲面拟合 2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下： ????? ??=k i i k x x x x B 0000 ， ????? ??=k j j k 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 T T ）（ 通过求解如下线性方程，即可得到系数矩阵C 。 UG B G G C B B T T T =)()( 2.2计算),(),,(* ***j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值 ),(**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)。 二、 源程序 ********* 第三次数值分析大作业Q9************ integer::i, j, K1, L1, n, m dimension X(31), Y(21), T(6), U(6), Z(6, 6), UX(11, 21), TY(11, 21), FXY(11, 21), C(6, 6) dimension z1(31, 21), z2(31, 21), z3(31, 21), z4(31, 21), z5(31, 21) dimension X1(8), Y1(5), FXY1(8, 5), PXY1(8, 5), UX1(8, 5), TY1(8, 5) 数值分析第二题 梁进明 SY0906529 算法设计方案。 一．矩阵的QR 分解。把矩阵A 分解为一个正交矩阵Q 与一个上三角矩阵R 的乘积，称为矩阵A 的 正三角分解，简称QR 分解。QR 分解的算法如下： 记1A A =，并记[]r ij n n Ar a ?=，令1Q I =（n 阶单位矩阵） 对于r=1,2,…,n-1执行 (1) 若(1,2,...,)r ir a i r r n =++全为零，则令1r r Q Q +=，1r r A A +=转(5);否则转(2) (2) 计算 2r d = ()sgn()r r rr r c a d =-(若() 0r rr a =,则取r r c d =) 2() r r r r rr h c c a =- (3) 令()()()1,(0,...,0,,,...,)r r r T n r rr r r r nr u a c a a R +=-∈ (4) 计算 11//r r r T r r r r r T r r r r T r r r r Q u Q Q u h p A u h A A u p ωω++==-==- (5) 继续 当此算法执行完后就得到正交矩阵n Q Q =和上三角矩阵n R A =且有A QR =。 二．矩阵的 拟上三角化。对实矩阵A 的拟上三角化具体算法如下： 记(1) A A =，并记()r A 的第r 列到第n 列的元素为(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于1,2,...,2r n =-执行 （1） 若() (2,3,...,)r ir a i r r n =++全为零，则令(1) ()r r A A +=,转(5);否则转(2)。 （2） 计算北航数值分析大作业一
数值分析作业答案
北航数值分析大作业第一题幂法与反幂法
北航数值分析大作业第二题精解
北航数值分析报告第三次大作业
上海大学_王培康_数值分析大作业
数值分析第二次大作业
北航数值分析大作业第二题
数值分析大作业QR分解
北航数值分析第二次大作业--QR分解
北航数值分析报告大作业第八题
数值分析大作业 三、四、五、六、七
北航数值分析课程第一次大作业讲解
北理工数值分析大作业
数值分析作业答案.doc
北航数值分析大作业3
北航数值分析大作业第二次
BUAA数值分析大作业三
北航数值分析大作业 第二题 QR分解
相关主题
文本预览