数值分析第二次作业
数值分析第二次作业及答案
数值分析第二次作业及答案1. 用矩阵的直接三角分解法（LU 分解）解方程组1234102050101312431701037x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦。
212223243132333441424344212223243132333441424344110201020101011124310103110201101,112110101l uu u l l u u l l l u luu u l l u u l l l u ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⇒=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎢⎥⎣⎦⎣⎦⎣⎦解：设111122223333444410201012125115102053101310136212117216420101724y y x x y y x x y y x x y y x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎥⎣⎦==⎡⎤⎧⎡⎤⎧⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎪⎢⎥⎪⎢⎥⎢⎥⎢⎥⎢⎥==⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⇒=⇒⎨⎨⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⎪⎪⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎪⎪==⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦⎣⎦⎩⎣⎦⎩由2.矩阵第一行乘以一数，成为211A λλ⎡⎤=⎢⎥⎣⎦,证明当23λ=±时，()cond A ∞有最小值。
112131213 21121223A AA A λλλλλλλλ--∞∞⎧⎛⎫≥- ⎪⎪⎛⎫⎪=⇒==⇒=+ ⎪⎨⎪ ⎪⎝⎭⎪<- ⎪⎪⎩⎝⎭解：123632() min ()7.22343cond A AA cond A λλλλλ-∞∞∞∞⎧+≥⎪⎪∴====⎨⎪+<⎪⎩故当时， 3. 设有方程组AX b =，其中121011221,,302223A b ⎡⎤⎢⎥-⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥-⎢⎥⎣⎦已知它有解12130X ⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦。
数值分析第二次大作业——编写ln函数实验报告
为一个整体，将 ln(a)进行 Taylor 级数分解，进而按式计算可以求得 ln(x)的值。由于我们选 取级数的前 n 项和近似 ln(a)，则有：
Rn (a
1)
(1)n
(a 1)n1 n 1
…；
|
Rn
(a
1)
||
(a 1)n1 n 1
|
。
2
自 92 乔晖 2009011414
精度为 20 位，即 h4 1020 。
int * GetNumber() {return number;}
//获得大数中的整数数组
bool GetSgn() const {return sgn;}
//获得大数的符号
bool IsZero();
//判断大数是否为 0
CBigNumber Reverse() {sgn = !sgn;return *this;} //取相反数
数值分析——编写 ln 函数
实验报告
自 92 乔晖 2009011414
自 92 乔晖 2009011414
数值分析——编写 ln 函数实验报告
目录
一、 实验任务 ................................................... 2
二、 编译环境 ................................................... 2
1u
x 1
35
5
自 92 乔晖 2009011414
数值分析——编写 ln 函数实验报告
其前 n 项和，近似可得 ln(x)的值。
基于以上原理，可以设计迭代算法，累加 2u2i1 得到前 n 项和。每次计算后，判断计 2i 1
数值分析第二次大作业
《数值分析》计算实习报告第二题院系：机械工程及自动化学院学号：姓名：2017年11月一、题目要求试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知a ij ={sin (0.5i +0.2j ) i ≠j1.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 ira i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。
（2）计算r d =()()1,sgn r r r r r c a d +=- (若()1,0r r r a +=，则取r r c d =)2()1,r r r r r r h c c a +=-（3）令()()()()1,2,0,,0,,,,T r r r n r r rr r rnru ac aaR ++=-∈。
数值分析大作业
第二次计算实验：SVD及其应用梁杰存2014310739航博1431.方法求矩阵A奇异值分解一个途径是求解A T A的特征值，但因为舍入误差容易丢掉小奇异值。
因此通常先将矩阵上双对角化，即构造正交阵Q和W,使得Q T AW=B(upper−bidiagnal)。
这一过程可以通过逐次Household变换或逐次Given’s变换完成，还有一种基于待定系数法思想的Lanczos算法。
由于Linpack中SVD算法需要输入上双对角矩阵，本文采用Lanczos 算法实现上双对角化。
1.1.隐式零移位QR法（implicit zero-shift QR）与传统移位QR迭代算法不同，隐式零移位QR算法不进行移位，并且第一步构造右乘Given’s变换矩阵GR(1,2)将上双对角矩阵B的(1,2)位置上的元素12b消零，而不是传统方法中引入一个非零元素21b。
但这一步可能会使原来为零的b12变为非零。
第二步左乘Given’s阵GL(1,2)使得12b为0，但可能会使为零b13变为非零。
与上述步骤类似，将b13变为0后可能会使b23非零。
如下图所示，重复上述步骤最终将恢复为上双对角矩阵，即完成一步隐式零移位QR迭代。
反复迭代，矩阵B将趋近与对角阵阵，对角元即特征值。
图1隐式QR迭代1.2.分而治之(Divide-and-conquer)分而治之算法将上双对角阵B分成有两个互相独立对角块矩阵与另一矩阵之和，即：B=B100B2+b m vvT=Q1Σ1Q1T00Q2Σ1Q2T+b m vv T =Q100Q2(Σ100Σ1+b m uuT)Q1T00Q2T所以矩阵B的特征值与矩阵D+ρu u T的特征值相同，其中D=Σ100Σ1为对角阵，又：det D+ρu u T−λI=det((D−λI)(I+ρD−λ−1u u T))由于D−λI非奇异，则det I+ρD−λ−1u u T=1+ρu T D−λ−1u=1+ρu i2d i−λ=0ni=1在每个d i与d i+1之间分布着一个特征值，可用牛顿法快速找到该特征值。
数值分析第二章答案
1．当1,1,2x =-时，()0,3,4f x =-,求()f x 的二次插值多项式。
解：0120121200102021101201220211,1,2,()0,()3,()4;()()1()(1)(2)()()2()()1()(1)(2)()()6()()1()(1)(1)()()3x x x f x f x f x x x x x l x x x x x x x x x x x l x x x x x x x x x x x l x x x x x x x ==-===-=--==-+-----==------==-+-- 则二次拉格朗日插值多项式为220()()k k k L x y l x ==∑0223()4()14(1)(2)(1)(1)23537623l x l x x x x x x x =-+=---+-+=+- 5设[]2(),f x Ca b ∈且()()0,f a f b ==求证： 21m ax ()()m ax ().8a x b a x bf x b a f x ≤≤≤≤''≤- 解：令01,x a x b ==，以此为插值节点，则线性插值多项式为10101010()()()x x x x L x f x f x x x x x --=+-- =()()x bx af a f b a b x a --=+--1()()0()0f a f b L x ==∴= 又 插值余项为1011()()()()()()2R x f x L x f x x x x x ''=-=--011()()()()2f x f x x x x x ''∴=--[]012012102()()1()()21()41()4x x x x x x x x x x b a --⎧⎫≤-+-⎨⎬⎩⎭=-=- 又 ∴21m ax ()()m ax ().8a x b a x bf x b a f x ≤≤≤≤''≤- 16．求一个次数不高于4次的多项式P （x ），使它满足(0)(0)0,(1)(1)0,(2)0P P P P P ''=====解：利用埃米尔特插值可得到次数不高于4的多项式0101010,10,10,1x x y y m m ======11300201001012()()()()(12)()(12)(1)j j j j j j H x y x m x x x x xx x x x x x x αβα===+--=---=+-∑∑210110102()(12)()(32)x x x x x x x x x x x α--=---=-2021()(1)()(1)x x x x x xββ=-=-22323()(32)(1)2H x x x x x x x ∴=-+-=-+设22301()()()()P x H x A x x x x =+--其中，A 为待定常数3222(2)1()2(1)P P x x x Ax x =∴=-++-14A ∴= 从而221()(3)4P x x x =-19．求4()f x x =在[,]a b 上分段埃尔米特插值，并估计误差。
研究生《数值分析》课程作业(二) (含答案)
研究生《数值分析》课程作业（二）姓名： 学号： 专业： 1、据如下函数值表，建立二次的Lagrange 插值多项式及Newton 插值多项式。
20012222()()()()()()()(1)(2)(0)(2)(-0)(1)593143(01)(02)(10)(12(20)(21)22L x f x l x f x l x f x l x x x x x x x x x =++-----=⨯+⨯+⨯=-+------解： 二次 l agr ange插值）Newton 插值多项式：200100120122()()[,](-)[,,](-)(-)555932(0)(0)(1)32()32222N x f x f x x x x f x x x x x x x x x x x x x x x =++=-⨯-+--=-+-=-+ ()y f x =2、已知单调连续函数在如下采样点处的函数值*()0[2,4],f x x =求方程在内根的近似值使误差尽可能小。
解：1()()y f x x f y -==解：对的反函数进行二次插值1110201122012010210122021(0)(0)(0)(0)(0)(0)(0)()()()()()()()()()(0 2.25)(05)(03)(05)(03)(0 2.25)2 3.54(3 2.25)(35)(2.253)(2.255)(53)(5 2.25)y y y y y y L f y f y f y y y y y y y y y y y y y ---------=++--------+-+-=⨯+⨯+⨯----+-+- 2.945≈()(1)01(1)1()[,]()(,),()[,],()()()()()(1)!,n n n n n n n n f x a b f x a b a x x x b L x x a b f R x f x L x x n a b x ξωξ+++≤<<<≤∈=-=+∈ 3、证明：设在上连续，在内存在，节点是满足拉格朗日插值条件的多项式，则对任何插值余项这里（）且依赖于。
数值分析第二次作业答案
练习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。
《应用数值分析》第二次上机作业
可利用递推关系
P ( x ) 1, P ( x ) x
0
1
P( n
x
) ( 2 n1)xP ( n1
x)
(n
1)P ( n2
x ) /
n
n 2, 3, .....
2、被逼近函数 f(x)不用内联函数构造，而改用 M 文件建立数学函数。这样，
此程序可通过修改建立数学函数的 M 文件以适用不同的被逼近函数（要
fa=feval(func,a);fb=feval(func,b); T2n=(fa+fb)*(b-a)/2.0; else n=2^(k-2); h=(b-a)/n; x=a+h/2; sum=0.0; for i=1:n
fx=feval(func,x); sum=sum+fx;
x=x+h;
end
二、被逼近函数用 M 文件建立（例如下面） function f=fun(x) f=1./(1+x.^2);
三、变步长复化梯形求积公式的算法
ba
1.h b a ,T
( f (a ) f (b ))
2
h 2.H 0, x a
2
3.H H f ( x ), x x h
2、三次最佳平方逼近函数：
输入 syms x;
[c,s]=leastp_ww(x,6,3,-1,1)
运行结果：
s=
17336850270911903/18014398509481984-9565132890499617/18014398509481984
*x^2
c=
0.785398004451730 0 -0.353981025621068
数值分析第二次大作业SOR最优松弛因子选取方法及SOR迭代法的改进
《数值分析》第二次大作业题目：SOR最优松弛因子选取方法及SOR迭代法的改进内容：1.SOR最优松弛因子选取方法2.SOR迭代法的改进（SSOR迭代法）3.SSOR迭代法的Matlab程序4.举例比较jacobi，Gauss-Seidel，SOR及SSOR 迭代法的收敛速度姓名：合肥工业大学学号：2011班级：信息与计算科学11-1班参考资料：1.《确定SOR最优松弛因子的一个实用算法》李春光等《计算力学学报》2.《数值分析与实验》，薛毅,北京工业大学出版社.3.《数值分析中的迭代法解线性方程组》，马云,科学出版社4.《非线性互补问题的改进超松弛迭代算法》，段班祥等,江西师范大学出版社5.《迭代法解线性方程组的收敛性比较》，郑亚敏,江西科学出版社.一、SOR最优松弛因子选取方法SOR迭代法迭代公式：x(k+1)i=(1-ω)xi+(k) bi-∑aijxjaii⎝j=1ω⎛ i-1(k+1)-j=i+1∑axijn(k)j⎫⎪ (i=1,2,..n.), ⎪⎭1.二分比较法将松弛因子1/2，ω的区间（1,2）进行二分，每个小区间的长度为ω去中间值3/2，按照SOR 迭代法迭代公式，求出跌代次数k,如果k不超过指定的发散常数，则可确定ω的值；否则将（1,2）四等分，每个区间长度为1/4,ω取各分点值，继续迭代，一般地，将1区间（1,2）二分M次，每次二分步长为，ω一次取取各分点值，2M按照SOR迭代法迭代公式，求出跌代次数k,如果k不超过指定的发散常数，则可确定的ω的值，这样总能找到一个不超过指定发散常数ω值。
2.逐步搜索法将1+ω的取值区间（1,2）进行M等分，ω分别取ω的值。
12M-1,1+,...,1+，通过迭代公式依次对同意精度要求求出迭代MMM次数k的值，并从中选出最优松弛因子3.黄金分割法依据黄金分割比的思想，通过计算机主动选取最优松弛因子的近似值，步骤如下a.对（1，2）区间进行第一次0.618的分割，区间边界a1=1,b1=2，在区间(a1,b1)分割出黄金点p1=a1+0.618(b1-a1)，进行SOR迭代法的迭代，求出迭代次数k的值，如果没有超过规定的发散常数，迭代结束，否则做步骤b。
数值分析报告第二题 北航 大作业
实用文档数值分析第二题目录数值分析第二题 (1)1 引言： (1)1.1 矩阵的拟上三角化 (1)1.2 矩阵的特征值求解 (2)1.3 矩阵的特征向量求解 (4)2 算法的程序实现 (5)2.1 主程序 (5)2.2 子程序的实现 (7)3计算结果 (8)3.2矩阵Q、R 以及乘积RQ (9)3.3 各实特征值及其相对应的特征向量 (10)4 实验结论 (12)附录源代码 (12)1 引言1.1 矩阵的拟上三角化为了减少计算量，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵 。
对A 拟上三角化，得到拟上三角矩阵)1(-n A ，具体算法如下： 记A A =)1(，并记)(r A 的第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)( ++=全为零，则令)()1(r r A A =+,转5；否则转2。
2. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若)(,12r rr r r r a c c h +-=3. 令()n Tr nrr 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 Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。
1.2 矩阵的特征值求解使用带双步位移的QR 方法计算矩阵)1(-n A 的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。
北航 数值分析第二次大作业(带双步位移的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)，再求解方程组。
数值分析第二次作业答案answer2
S4 = 0.11157238253891，S8 = 0.11157181325263。 同学们根据自己理解计算 S4 ，S8 都可。 复合梯形公式和复合 Simpson 公式的代码已重复多次，同学们自己整 理。 3. 用 Simpson 公式计算积分 误 差 为 |R(f )| = | − η ∈ (0, 1)。 4. 推导下列三种矩形求积公式： ∫b f (x)dx ∫a b f (x)dx ∫a b a f (x)dx = (b − a)f (a) + = (b − = (b −
14.7 53.63 从而 a = −7.855048，b = 22.25376。 2. 已知实验数据如下： 。 xi 19 25 31 38
44
yi 19.0 32.3 49.0 73.3 97.8 用最小二乘法求形如 y = a + bx2 的经验公式。 答案：两个待定常数，只能两个 φ。 φ0 ，φ1 也必须形如 y = a + bx2 。 可设 φ0 = 1，φ1 = x2 。法方程为： ( 5 5327 )( a b ) = ( 271.4 369321.5 )
第三章 函数逼近 1. 观测物体的直线运动，得出以下数据： 时间 t(s) 0 0.9 1.9 3.0 3.9 5.0 距离 s(m) 0 求运动方程。 ( 10 φ0 = 1，φ1 = t。法方程为： 6 14.7 )( a b ) = ( 280 1078 )
6
1. 用 LU 分解及列主元高斯消去法解线性方程组 8 10 −7 0 1 x1 −3 2.099999 6 2 x 5.900001 2 = 5 5 − 1 5 − 1 x 3 x4 1 2 1 0 2 输出 Ax = b 中系数 A = LU 分解的矩阵 L 及 U ，解向量 x 及 det A；列 主元法的行交换次序，解向量 x 及 det A；比较两种方法所得的结果。 代码： A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2]; b=[8,5.900001,5,1]'; x=A\b;x(1) 结果：1.7764e-016 LU分解代码： A=[10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2]; b=[8,5.900001,5,1]'; [m,n] = size(A); if m~=n, error('A matrix needs to be square'); end for i=1:n-1 pivot = A(i,i); if abs(pivot)<50*eps, error('zero pivot encountered'); end for k = i+1:n A(k,i) = A(k,i)/pivot; A(k,i+1:n) = A(k,i+1:n) - A(k,i)*A(i,i+1:n); end end 7
高等数值分析第二次实验作业(贾忠孝)
3. 对 1 中的例子固定特征值的实部, 变化虚部, 比较收敛性. A 矩阵每个小块有如下形式：
A i i
其特征值为
i i
i ii ，对其次对角元素 i 做如下处理：
i ' k i
其特征值变为 果：
i ik i 。在本题中分别设定 k 的取值为 0.5，1，2，4。可以得到下面的结
结果分析： （1） 当 k 越大，既即特征根虚部越大时，基本 Arnoldi 方法和 GMRES 方法的收敛速度会 越来越慢，Arnoldi 方法的振荡程度会随着 k 值的增大而增大，振荡范围也增大。 当 k 大到 10 或者 10 以上时，两法都只在最后一步（第 1000 步）收敛，而且都是 急剧收敛。 （2） GMRES 法的残量一直比 Arnoldi 法的残量小，且相对残差的曲线平滑。.
GMRES算 法 计 算 m个 特 征 向 量 组 成 的 b的 收 敛 曲 线 4 2 0 -2 GMRES, GMRES, GMRES, GMRES, m m m m = = = = 10 50 100 1000
log(||rk||)
-4 -6 -8 -10 -12 -14
0
50 Iteration Times
100
150
结果分析： （1） 当 m 较小时，收敛速度很快，但其随着 m 值的增大而慢慢变小。但是一直到最后 一步之前，收敛的幅度都不大，往往是在最后一步急剧收敛； （2） 当 m 较大时，收敛的速度随 m 的增大而微弱减小，而且其在步数比较小的时候， 收敛速度较快，越往后走反而慢。 （3） Arnoldi 法比 GMRES 法收敛曲线的趋势是基本一样的，而且 Arnoldi 法收敛比 GMRES 法要快，相对残量： 1.构造例子特征值全部在右半平面时,观察基本的 Arnoldi 方法和 GMRES 方法的数值性态,和 相应重新启动算法的收敛性. 解：构造 1000 阶符合条件的矩阵 A。根据实 Schur 分解，构造如下形式的矩阵：
数值分析作业(第二章)
数值分析作业（第二章）习题1. 当x=1,-1,2时，f(x)=0,-3,4,求()f x 的二次插值多项式。
（1）用单项式基底； （2）用拉格朗日插值基底； （3）用牛顿基底。
证明三种方法得到的多项式是相同的。
解：（1）假设f(x)的二次插值多项式为：2210()f x a x a x a =++由于 x=1,-1,2时，f(x)=0,-3,4则有 2100a a a ++=；2103a a a -+=-；210424a a a ++=求得 256a =；132a =；073a =- 则有 2537()623f x x x =+-（2）用拉格朗日插值基底： 由于 0121,1,2,x x x ==-=012()0,()3,()4;f x f x f x ==-=则有 1200102()()1()(1)(2)()()2x x x x l x x x x x x x --==-+---0211012()()1()(1)(2)()()6x x x x l x x x x x x x --==----0122021()()1()(1)(1)()()3x x x x l x x x x x x x --==-+--拉格朗日插值多项式为：220011220()()()()()()()()k k k L x y l x f x l x f x l x f x l x ===++∑2537623x x =+- 则f(x)二次插值多项式为：22537()623L x x x =+-（3）采用牛顿基底： 均差表如下所示：则有牛顿插值多项式为：[][]2001001201()(),(),,()()N x f x f x x x x f x x x x x x x =+-+-- 2537623x x =+- 则f(x)二次插值多项式为：22537()623L x x x =+- 由以上计算可知，三种方法得到的多项式是相同的。
数值分析
数值分析实验课第二次作业数学与统计学院12级师范二班1207410056 李军旗第一题1.编写命令文件：输入向量a和向量b，将a和b分别进行从小到大排序并赋值给c和d，计算并输出c+d的值。
>> a=input('a=');b=input('b=');(输入)c=sort(a);d=sort(b);（向量大小排序）s=c+d;a=[4,2,1];b=[5,3,7];（任意输入向量，但必须保证a和b维数相同）>> ss =4 7 11第二题2.求Fibonacci数列(1)大于4000的最小项。
(2)5000之内的项数。
function f=f(n)f=[1,1];i=1;while f(i)+f(i+1)<nf(i+2)=f(i)+f(i+1);i=i+1;end (定义函数并保存)>>f(5000) （调用函数并输出）ans =Columns 1 through 81 123 5 8 13 21Columns 9 through 1634 55 89 144 233 377 610 987Columns 17 through 191597 2584 4181（由上可知大于4000的最小项是4181,5000内一共19项）第三题3.求给定矩阵A中的所有素数的和。
>> A=input('输入矩阵A=')[m,n]=size(A);prime=[]; （保存矩阵的素数）sum=0;B=0;for k=1:mfor i=1:nif A(k,i)==2|A(k,i)==3;B=1; （2,3素数单独拿出）endfor j=2:fix(sqrt(A(k,i)))if mod(A(k,i),j)==0B=0;break; （如果是素数记为1，不是记为0）elseB=1endendif B==1prime=[prime,A(k,i)]; （选出记为1的素数求和） sum=sum+A(k,i);endendendprimesum输入矩阵A=[1,3,5;2,4,6;7,9,11]A =1 3 52 4 67 9 11B =1B =1B =1B =1B =1prime =3 5 2 7 11sum =28第四题4.编写函数文件piecef.m, 通过调用piecef（x）计算函数在x点的值，该函数是如下形式的分段函数，绘制出其在[0，10]上的曲线。
数值分析-第二次作业答案
2020级数值分析第二次作业(非线性方程求根)参考答案和评分标准班级学号姓名一(20分)用二分法求方程3()10f x x x =--=在区间[1.0,1.5]内的一个实根，且要求有3位有效数字。
试完成：(1)估计需要二分的次数；(8分)解：容易知道方程在[1.0,1.5]有且仅有一个实根。
记此实根为*x ，根据二分法误差估计公式有12()1*22k k k b a x x ++--≤=要使得近似解有3位有效数字，只需要有22111022k -+≤⨯从而可得6k ≥，即满足精度要求的二分次数为6次。
(2)将计算过程中数据填入表1.(中间过程填写到小数点后面3位)(12分，每个k x 得2分，其它空不计分)表1题1计算过程kk a k b kx 0 1.0 1.5 1.251 1.25 1.5 1.3752 1.25 1.375 1.3163 1.313 1.375 1.3494 1.313 1.344 1.3285 1.313 1.328 1.32061.3201.3281.324二.(10分)为了计算方程()3sin 2120f x x x =--=的根，某同学将()0f x =改写为14sin 23x x =+，并建立迭代公式114sin 23k k x x +=+。
请问此迭代公式在R 上是否全局收敛的吗？说明理由。
证明：(1)对任意的x R ∈，有11113()4sin 2,333x x R ϕ⎡⎤=+∈⊆⎢⎣⎦；(2)对任意的x R ∈，有22'()cos 2133x x ϕ=≤<;从而可知，迭代格式在R 上全局收敛。
三.(20分)设有方程3()10f x x x =--=，试回答下列问题：(1)确定方程3()10f x x x =--=实根的数目；(4分)解：由2'()31f x x =-可知函数3()1f x x x =--的单调递增区间是,33⎛⎛⎫-∞-⋃+∞ ⎪ ⎪⎝⎭⎝⎭，单调递减区间是,33⎛- ⎝⎭。
北航数值分析大作业第二题
北航数值分析大作业第二题本页仅作为文档封面，使用时可以删除This document is for reference only-rar21year.March数值分析第二次大作业史立峰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. 计算()∑+==nr i r irr a d 12)(()()r r r r r r r r r r d c a d a c ==-=++则取,0sgn )(,1)(,1若 )(,12r rr r r r a c c h +-=3. 令()nTr nrr 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 Tr r h u p t /=r r r r u t q -=ωT rr T r r r r p u u A A --=+ω)()1(5. 继续。
（3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下：1. 给定精度水平0>ε和迭代最大次数L 。
西南交大数值分析第二次大作业(可以运行)
数值分析第二次大作业1（1）用Lagrange插值法程序：function f=Lang(x,y,x0)syms t;f=0;n=length(x);for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));endf=f+l;simplify(f);if(i==n)if(nargin==3)f=subs(f,'t',x0);elsef=collect(f);f=vpa(f,6);endendendx=[1,2,3,-4,5];y=[2,48,272,1182,2262]; t=[-1];disp('插值结果')yt=Lang(x,y,t)disp('插值多项式')yt=Lang(x,y)ezplot(yt,[-1,5]);运行结果：插值结果：Yt= 12.0000插值多项式：yt =4.0*t^4 - 2.0*t^3 + t^2 - 3.0*t + 2.0（2）构造arctan x在[1,6]基于等距节点的10次插值多项式程序：function f=New(x,y,x0)syms t;if(length(x)==length(y))n=length(x);c(1:n)=0.0;elsedisp('xºÍyάÊý²»µÈ£¡');return;endf=y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1)));for(i=1:n-1)for(j=1:n-i)y1(j)=y(j+1)-y(j);endc(i)=y1(1);l=t;for(k=1:i-1)l=l*(t-k);end;f=f+c(i)*l/factorial(i);simplify(f);y=y1;if(i==n-1)if(nargin==3)f=subs(f,'t',(x0-x(1))/(x(2)-x(1)));elsef=collect(f);f=vpa(f,6);endendend>>x=[1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6];y=[atan(1),atan(1.5),atan(2),atan(2.5),atan(3),atan(3.5),atan(4),atan(4.5),atan(5),atan (5.5),atan(6)];disp('插值多项式')yt=New(x,y)ezplot(yt,[1,6]);hold onezplot('atan(t)',[1,6])grid on运行结果：插值多项式yt =1.34684*10^(-10)*t^10 - 6.61748*10^(-9)*t^9 + 1.28344*10^(-7)*t^8 - 0.00000104758*t^7 - 0.00000243837*t^6 + 0.000149168*t^5 - 0.00176296*t^4 + 0.0125826*t^3 - 0.0640379*t^2 + 0.250468*t + 0.7853982（1）用MATLAB自带spline函数用于进行三次样条插值程序：>>x=[-5,-4,-3,-2,-1,0,1,2,3,4,5];y=[0.03846,0.05882,0.10000,0.20000,0.50000,1.00000,0.50000,0.20000,0.10000,0.0 5882,0.03846];xi=linspace(-5,5)yi=spline(x,y,xi);plot(x,y,'rp',xi,yi);hold on;syms xfx=1/(1+x^2);ezplot(fx);grid on运行结果：由图可知，三次样条插值多项式图像与原函数图像基本一致。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析大作业（二）院系：学号：姓名：一、算法设计方案：1、矩阵A的存储和特征值的数据结构设计A是一个10×10的方阵，且由给定的元素计算公式知A是非对称的矩阵，因此采取二维数组来存储数组A。
考虑到特征值为复数的情况，采用结构体来存储特征值的实部和虚部，当特征值为实数时，虚部赋值为零。
2、求解过程的分析首先利用HouseHolder矩阵将A拟上三角化，得到A(n-1)。
然后利用带双步位移的QR分解法求解矩阵A的特征值。
利用选列主元素的高斯消去法求解A的实特征值对应的特征向量，将等式两边移项后，相当于求解齐次线性方程组。
在计算前首先要求齐次方程组的系数矩阵的秩，从而确定自由量,赋予适当的初值，从而可以求得特征值对应的一个特征值。
利用QR分解的基本方法可以得到A(n-1)进行QR分解后的Q、R和RQ。
二、源程序#include <stdio.h>#include <math.h>#define N 10#define epsilon 1E-12#define L 10000typedef struct Comp{double real;double image;}Comp;void getRealMatrix(int n,double A[][n],double Q[][n],double regulate); int isContainThisValue(int m,int mainIndex[],int index);int maxElementIndex(int m,int n,int colNum,double array[][n],int mainIndex[]);void mainElementElimination(int n,double equationArray[][n],double solveArray[],double);void qrMethod(int,double a[][]);int sgn(double realNum);void outputMatrix(int m,int n,double A[][n]);void calcMatrix(int m,int n,double A[][n],double M[][n],double s,double t);void productMatrix(int m,int n,double A[][n],double M[][n]);void hessenbergMatrix(int n,double matrix_a[][]);void solvePowerEquation(double b,double c,Comp twoDArr[2]);void qrIterateFormula(int m,int n,double M[][n],double A[][n]);void initialMatrix(int n, double A[][]);double crossProduct(double *p_1,double *p_2,int startIndex,int endIndex,int increment);void transpositionMatirx(int m,int n,double matrix[][]);void iterateFormula(int m,int n,double matrix[][],double w[],double u[],double p[]);void qrWithTwoStepMove(int n,double matrix[][],double matrix_Q[][],Comp lamdaValue[],double e,double iterCount);int main(){double Matrix_Q[N][N],Matrix_A[N][N];Comp lamdas[N];double vector[N]={0};int i,j;//初始化矩阵initialMatrix(N,Matrix_A);//拟上三角化化hessenbergMatrix(N,Matrix_A);//输出拟上三角化后的矩阵printf("A拟上三角化的矩阵：\n");outputMatrix(N,N,Matrix_A);//双步位移的QRqrWithTwoStepMove(N,Matrix_A,Matrix_Q,lamdas,epsilon,L);//求解Q、R、RQinitialMatrix(N,Matrix_A);hessenbergMatrix(N,Matrix_A);qrMethod(N,Matrix_A);for(i=0;i<=N-1;i++){printf("%.12E,%.12E\n",lamdas[i].real,lamdas[i].image);}//求解实特征值对应的特征向量initialMatrix(N,Matrix_A);for(i=0;i<=N-1;i++){if(lamdas[i].image!=0)continue;else{getRealMatrix(N,Matrix_A,Matrix_Q,lamdas[i].real);mainElementElimination(N,Matrix_Q,vector,1);printf("%.12E的特征向量为\n",lamdas[i].real);for(j=0;j<=N-1;j++)printf("%.12E ",vector[j]);printf("\n");}}return 0;}//输出矩阵void outputMatrix(int m,int n,double A[][n]){int i,j;for(i=0;i<=m-1;i++){for(j=0;j<=m-1;j++){printf("%0.12E ",A[i][j]);}printf("\n");}printf("\n");}//初始化矩阵的元素void initialMatrix(int n, double A[][n]){int i,j;for(i=0;i<=n-1;i++)for(j=0;j<=n-1;j++){if(i==j)A[i][j] = 1.5*cos((i+1)+1.2*(j+1));elseA[i][j] = sin(0.5*(i+1)+0.2*(j+1));}}//对矩阵进行拟上三角化void hessenbergMatrix(int n,double matrix_a[][n]){int r,i;double u[n],p[n],w[n],c,d,h;for(r=0;r<=n-3;r++){if(crossProduct(&matrix_a[r+2][r],&matrix_a[r+2][r],r+2,n-1,n) == 0)continue;d = sqrt(crossProduct(&matrix_a[r+1][r],&matrix_a[r+1][r],r+1,n-1,n)); c = matrix_a[r+1][r] == 0?d:(-sgn(matrix_a[r+1][r])*d);h = c*c - c*matrix_a[r+1][r];for(i=0;i<=n-1;i++){u[i] = (i<(r+1))?0:(i==(r+1))?(matrix_a[r+1][r]-c):matrix_a[i][r];}transpositionMatirx(n,n,matrix_a);for(i=0;i<=n-1;i++){p[i] = crossProduct(&matrix_a[0][0]+n*i,u,0,n-1,1)/h;}transpositionMatirx(n,n,matrix_a);for(i=0;i<=n-1;i++){w[i] = (crossProduct(&matrix_a[0][0]+n*i,u,0,n-1,1)- crossProduct(p,u,0,n-1,1)*u[i])/h;}iterateFormula(n,n,matrix_a,w,u,p);}}//双步位移的QR算法实现void qrWithTwoStepMove(int n,double matrix[][n],double matrix_Q[][n],Comp lamdaValue[],double e,double iterCount){int k,m,i=0,index=0;double s,t;Comp solver[2];k =1;m = n-1;while(1){if(fabs(matrix[m][m-1])<=e){lamdaValue[index].real = matrix[m][m];lamdaValue[index].image =0;index++;m=m-1;if(m<2){if(m==0){lamdaValue[index].real= matrix[m][m];lamdaValue[index].image = 0;index++;}else if(m==1){s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];solvePowerEquation(s,t,solver);for(i=0;i<=1;i++,index++)lamdaValue[index] = solver[i];}return;}}else{if(fabs(matrix[m-1][m-2])<=e){s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];solvePowerEquation(s,t,solver);for(i=0;i<=1;i++,index++)lamdaValue[index] = solver[i];m=m-2;if(m<2){if(m==0){lamdaValue[index].real= matrix[m][m];lamdaValue[index].image = 0;index++;}else if(m==1){s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];solvePowerEquation(s,t,solver);for(i=0;i<=1;i++,index++)lamdaValue[index] = solver[i];}return;}}else{if(k == iterCount)return;else{s = matrix[m-1][m-1]+ matrix[m][m];t = matrix[m-1][m-1]*matrix[m][m] - matrix[m][m-1]*matrix[m-1][m];calcMatrix(m,n,matrix,matrix_Q,s,t);qrIterateFormula(m+1,n,matrix_Q,matrix);k =k + 1;}}}}}//矩阵拟上三角化的迭代式子void iterateFormula(int m,int n,double matrix[][n],double w[],double u[],double p[]){int i,j;for(i=0;i<=m-1;i++)for(j=0;j<=m-1;j++){matrix[i][j] = matrix[i][j]-w[i]*u[j] - u[i]*p[j];}}//计算向量的内积， startIndex是起始下标，endIndex是结束下标，increment是元素增量double crossProduct(double *p_1,double *p_2,int startIndex,int endIndex,int increment){double sum = 0;int i;if(p_1 == p_2){for(i=startIndex;i<=endIndex;i++){sum += (*p_1)*(*p_1);p_1 = p_1 + increment;}}else{for(i=startIndex;i<=endIndex;i++){sum += (*p_1)*(*p_2);p_1 = p_1 + increment;p_2 = p_2 + increment;}}return sum;}//矩阵转置void transpositionMatirx(int m,int n,double matrix[][n]){int i,j;double changeVar;for(i=0;i<=m-1;i++){for(j=i;j<=m-1;j++){if(i != j){changeVar = matrix[i][j];matrix[i][j] = matrix[j][i];matrix[j][i] = changeVar;}}}}//符号函数int sgn(double realNum){return realNum== 0? 0 : (realNum>0?1:-1);}//代入的矩阵void calcMatrix(int m,int n,double A[][n],double M[][n],double s,double t){int i,j;productMatrix(m,n,A,M);for(i=0;i<=m;i++)for(j=0;j<=m;j++){M[i][j] = M[i][j] - s*A[i][j] + ((i==j)? t:0);}}//计算A平方void productMatrix(int m,int n,double A[][n],double M[][n]){int i,j,r;double sum;for(i=0;i<=m;i++)for(j=0;j<=m;j++){sum = 0;for(r=0;r<=m;r++)sum += A[i][r]*A[r][j];M[i][j]= sum;}}//计算lamdavoid solvePowerEquation(double b,double c,Comp twoDArr[2]){double delta = b*b-4*c;if(delta<0){twoDArr[0].real = b/2;twoDArr[0].image = sqrt(fabs(delta))/2;twoDArr[1].real = b/2;twoDArr[1].image = -sqrt(fabs(delta))/2;}else{twoDArr[0].real = (b+sqrt(delta))/2;twoDArr[0].image = 0;twoDArr[1].real = (b-sqrt(delta))/2;twoDArr[1].image = 0;}}//双步移位QR中的迭代公式void qrIterateFormula(int m,int n,double M[][n],double A[][n]) {int r,i;double c,d,h,u[m],v[m],p[m],w[m];for(r=0;r<=m-2;r++){if(crossProduct(&(M[r+1][r]),&(M[r+1][r]),r+1,m-1,n) == 0) {continue;}else{d = sqrt(crossProduct(&M[r][r],&M[r][r],r,m-1,n));c = (M[r][r] == 0)?d:(-sgn(M[r][r])*d);h = c*c - c*M[r][r];for(i=0;i<=m-1;i++){u[i] = (i<r)?0:(i==r)?(M[r][r]-c):M[i][r];}transpositionMatirx(m,n,M);for(i=0;i<=m-1;i++){v[i] = crossProduct(&M[0][0]+n*i,u,0,m-1,1)/h;}transpositionMatirx(m,n,M);for(i=0;i<=m-1;i++)p[i]=0;iterateFormula(m,n,M,u,v,p);//第二部分计算transpositionMatirx(m,n,A);for(i=0;i<=m-1;i++){p[i] = crossProduct(&A[0][0]+n*i,u,0,m-1,1)/h;}transpositionMatirx(m,n,A);for(i=0;i<=m-1;i++){w[i] = (crossProduct(&A[0][0]+n*i,u,0,m-1,1) - crossProduct(p,u,0,m-1,1)*u[i])/h;}iterateFormula(m,n,A,w,u,p);}}}//基本QR方法void qrMethod(int n,double A[][n]){int r,i,j;double c,d,h,sum;double u[n],w[n],p[n];double Q[n][n];double R[n][n];//初始化为单位阵for(i=0;i<=n-1;i++)for(j=0;j<=n-1;j++){if(i==j){Q[i][j]=1;R[i][j] = 0;}else{Q[i][j]=0;R[i][j] = 0;}}for(r=0;r<=n-1;r++){//判断第r列元素从 r+1行开始是否全部为0if(crossProduct(&A[r+1][r],&A[r+1][r],r+1,n-1,n)==0)continue;else{d = sqrt(crossProduct(&A[r][r],&A[r][r],r,n-1,n));c = (A[r][r] == 0)?d:(-sgn(A[r][r])*d);h = c*c - c*A[r][r];//初始化向量ufor(i=0;i<=n-1;i++){u[i] = (i<r)?0:(i==r)?(A[r][r]-c):A[i][r];}for(i=0;i<=n-1;i++){w[i] = crossProduct(&Q[0][0]+n*i,u,0,n-1,1)/h;}//辅助作用for(i=0;i<=n-1;i++)p[i]=0;//迭代出QiterateFormula(n,n,Q,w,u,p);transpositionMatirx(n,n,A);for(i=0;i<=n-1;i++){p[i] = crossProduct(&A[0][0]+n*i,u,0,n-1,1)/h;}//辅助作用for(i=0;i<=n-1;i++)w[i]=0;transpositionMatirx(n,n,A);iterateFormula(n,n,A,w,u,p);}}printf("正交阵Q: \n");outputMatrix(n,n,Q);printf("上三角阵R: \n");outputMatrix(n,n,A);for(i=0;i<=n-1;i++){for(j=0;j<=n-1;j++){sum=0;for(r=i;r<=n-1;r++)sum += A[i][r]*Q[r][j];R[i][j] = sum;}}printf("输出RQ：\n");outputMatrix(N,N,R);}//得到齐次线性方程组的系数矩阵void getRealMatrix(int n,double A[][n],double Q[][n],double regulate) {int i,j;for(i=0;i<=n-1;i++)for(j=0;j<=n-1;j++){if(i==j)Q[i][j] = A[i][j] - regulate;elseQ[i][j] = A[i][j];}}//列主元素消元法void mainElementElimination(int n,double equationArray[][n],double solveArray[],double referValue){int mainIndex[n];int temp;int k=0,i=0,j=0;double sum = 0;double ratio,afterDivide;for(i=0;i<=n-1;i++)mainIndex[i] = -1;for(k=0;k<=n-1;k++){//寻找最大元素所在的行号mainIndex[k]=maxElementIndex(n,n,k,equationArray,mainIndex);temp = mainIndex[k];for(i=0;i<=n-1;i++){if(isContainThisValue(n,mainIndex,i))continue;else{ratio = equationArray[i][k]/equationArray[mainIndex[k]][k];for(j=k;j<=n-1;j++){equationArray[i][j] = equationArray[i][j] - equationArray[mainIndex[k]][j]*ratio;afterDivide = equationArray[i][j];}}}}//代入一个设定值solveArray[n-1] = referValue;for(i=n-2;i>=0;i--){sum = 0;for(j=i+1;j<=n-1;j++){sum += equationArray[mainIndex[i]][j]*solveArray[j];}solveArray[i] = - sum/equationArray[mainIndex[i]][i];}}//返回一列中最大元素的下标int maxElementIndex(int m,int n,int colNum,double array[][n],int mainIndex[]){int maxIndex=0;int i=0,j=0;double max = array[0][colNum];for(i=0;i<=m-1;i++){if(!isContainThisValue(m,mainIndex,i)){max = array[i][colNum];maxIndex = i;break;}}for(i++;i<=m-1;i++){if(isContainThisValue(m,mainIndex,i))continue;if(max<array[i][colNum]){max = array[i][colNum];maxIndex = i;}}return maxIndex;}//判断index是否在mainIndex数组中int isContainThisValue(int m,int mainIndex[],int index) {int i = 0;for(i=0;i<=m-1;i++){if(index == mainIndex[i])return 1;}return 0;}三、计算结果。