数值分析第二次大作业
数值分析第二次作业及答案
数值分析第二次作业及答案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 ⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦。
北航数值分析大作业二(纯原创,高分版)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令 以及最大迭代步数L；
2.若m≤0，则结束计算，已求出A的全部特征值，判断 或 或m≤2是否成立，成立则转3，否则转4；
3.若 ，则得一个特征值 ，m=m-1，降阶；若 ，则计算矩阵：
的特征值得矩阵A的两个特征值，m=m-2，降阶，转2.；
4.若k≤L，成立则令
k=k+1，转2，否则结束计算，为计算出矩阵A的全部特征值；
数值分析第二次大作业——编写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
北航数值分析全部三次大作业
北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。
我们被要求实现各种常用的线性方程组求解算法，例如高斯消元法、LU分解法和迭代法等。
我首先学习了这些算法的原理和实现方法，并借助Python编程语言编写了这些算法的代码。
在实验中，我们使用了不同规模和条件的线性方程组进行测试，并比较了不同算法的性能和精度。
通过这个作业，我深入了解了线性方程组求解的原理和方法，提高了我的编程和数值计算能力。
第二次大作业是关于数值积分的方法。
数值积分是数值分析中的重要内容，它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。
在这个作业中，我们需要实现不同的数值积分算法，例如矩形法、梯形法和辛普森法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们计算了不同函数的积分值，并对比了不同算法的精度和效率。
通过这个作业，我深入了解了数值积分的原理和方法，提高了我的编程和数学建模能力。
第三次大作业是关于常微分方程的数值解法。
常微分方程是数值分析中的核心内容之一，它可以用于描述众多物理、化学和生物现象。
在这个作业中，我们需要实现不同的常微分方程求解算法，例如欧拉法、龙格-库塔法和Adams法等。
我学习了这些算法的原理和实现方法，并使用Python编写了它们的代码。
在实验中，我们解决了一些具体的常微分方程问题，并比较了不同算法的精度和效率。
通过这个作业，我深入了解了常微分方程的原理和方法，提高了我的编程和问题求解能力。
总的来说，北航数值分析课程的三次大作业非常有挑战性，但也非常有意义。
通过这些作业，我在数值计算和编程方面得到了很大的提升，也更加深入地了解了数值分析的理论和方法。
虽然这些作业需要大量的时间和精力，但我相信这些努力将会对我未来的学习和工作产生积极的影响。
研究生《数值分析》课程作业(二) (含答案)
研究生《数值分析》课程作业（二）姓名： 学号： 专业： 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、证明：设在上连续，在内存在，节点是满足拉格朗日插值条件的多项式，则对任何插值余项这里（）且依赖于。
数值分析大作业之2
数值分析大作业二一、算法设计方案(一)算法流程1.将要求解的矩阵进行Householder变换，进行拟上三角化，并输出拟上三角化的结果；2.将拟上三角化后的矩阵进行带双步位移的QR分解(最大迭代次数1000次)，逐个求出特征值，输出QR法结束后得到的Q、R、RQ阵，输出求出的特征值（用实部和虚部表示）3.对于求出来的实特征值，使用带原点平移的反幂法求出对应的特征向量，输出这些特征向量。
（二）程序设计流程1. 定义精度和最大迭代次数；2. 使用容器vector进行编程（方便增减元素），使用传引用调用（提高执行效率）；3. 将各个步骤用到的数学算法封装成函数，方便调用。
具体需要的函数如下：double VectMod(vector< double > &p)：求向量的模vector< double > NumbMultiVect(vector< double > &vect, double a)：数乘向量double VectMultiVect(vector< double > &y, vector< double > &u)：求两个向量的积vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector<double > &u)：矩阵的转置乘向量vector< double > ArraMultiVect(vector< vector< double > > &B, vector< double > &u)：矩阵乘向量vector< vector< double > > VectMultiCovVect(vector< double > &a, vector< double >&b)：向量乘向量转置得到矩阵void ArraSubs(vector< vector< double > > &C, vector< vector < double > > D) ：两个矩阵相减Vector< double > VectSubs(vector< double > &a, vector< double > &b)：两个向量相减vector<vector<double>> GetM(vector<vector<double>> &A)：求解矩阵M （用于双部位移QR迭代用）double GetMode(vector < vector < double > > &B, const int r)：求解矩阵B的r列向量的模double GetModeH(vector < vector < double > > &B, const int r)：求解矩阵B的第r列向量的模，用于拟对角化vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)：一个实数乘矩阵bool IsBirZeroH(vector< vector< double > > &B, const int r)：判断B[i][r]对角线下是否为零void GausElim(vector< vector< double > > a)：列主元高斯消元法求齐次方程解向量void Stop(vector< vector< double > > &Ar)：停止，结束程序void SolutS1S2(complex< double > &s1, complex< double > &s2, vector< vector<double > > &A)：求解二阶子阵的特征值s1，s2；void Save2(complex< double > &s1, complex< double > &s2)：保存两个特征值void Save1(complex< double > &s)：保存一个特征值void JudgemBelow2(vector< vector< double > > &A, vector< vector< double > > Abk)：对于m == 1 及m == 0 的处理void Hessenberg(vector< vector< double > > &A)：矩阵拟上三角化void QRMethod(vector< vector< double > > A)：矩阵QR分解void CalculatAk(vector< vector < double > > &Ak)：带双步位移QR迭代法二、源程序#include "stdafx.h"#include <vector>#include <iostream>#include <math.h>#include <complex>#include <fstream>using namespace std;const double epsion = 1e-12;const int L = 1000;int m,n;int k = 1;vector< vector< double > > I;vector< complex< double > > Lambda;///////////////////////////以下为自定义的算法流程中用到的函数double VectMod(vector< double > &p) //求向量的模{double value = 0.0;vector< double >::size_type i,j;j = p.size();for (i=1; i<j; i++){value += p[i] * p[i];}value = sqrt(value);return value;}vector< double > NumbMultiVect(vector< double > &vect, double a) //数乘向量{int j = vect.size();vector< double > b(j, 0);for (int i=1; i<j; i++){b[i] = a * vect[i];}return b;}double VectMultiVect(vector< double > &y, vector< double > &u)//两个向量相乘{vector< double >::size_type a = y.size();double value = 0;for (vector< double >::size_type i=1; i<a; i++){value += y[i] * u[i];}return value;}vector< double > ConveArraMultiVect(vector< vector< double > > &B, vector< double > &u)//矩阵的转置乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[j][i] * u[j];}}return vec;}}vector< double > ArraMultiVect(vector< vector< double > > &B,vector< double > &u) //矩阵乘向量{int a = B.size();int b = u.size();vector< double > vec(a, 0);if (a != b){cerr << "Array and Vector not match in size!";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){vec[i] += B[i][j] * u[j];}}return vec;}}vector< vector<double> > GetM(vector< vector <double> > &A){int a = A.size();double s = A[a-2][a-2] + A[a-1][a-1];double t = A[a-2][a-2] * A[a-1][a-1] - A[a-1][a-2] * A[a-2][a-1];vector<vector<double>> D(a, vector< double >(a, 0));for (int i=1; i<a; i++){{double sum = 0;for (int k=1; k<a; k++){sum += A[i][k] * A[k][j];}D[i][j] = sum - s * A[i][j] + t *(i==j ? 1.0 : 0);}}return D;}double GetMode(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}double GetModeH(vector < vector < double > > &B, const int r){double value = 0;int a = B.size();for (int k=r+1; k<a; k++){value += B[k][r] * B[k][r];}value = sqrt(value);return value;}vector< vector< double > > NumbMultiArra(vector< vector< double > > &D, double a)//数乘//向量{int b = D.size();vector< vector< double > > U(b, vector< double >(b, 0));for (int i=1; i<b; i++){{U[i][j] = a *D[i][j];}}return U;}bool IsBirZero(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+1; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}bool IsBirZeroH(vector< vector< double > > &B, const int r){bool b = true;int a = B.size();for (int i=r+2; i<a; i++){if(abs(B[i][r]) > epsion){b = false;}}return b;}vector< vector< double > > VectMultiCovVect(vector< double > &a,vector< double > &b)//向量乘向量的转置得到矩阵{int s1 = a.size();int s2 = b.size();if (s1 != s2){cerr << "Vectors not match in size ! ";}else{vector< vector< double > > U(s1, vector< double >(s1, 0));for (int i=1; i<s1; i++){for (int j=1; j<s1; j++){U[i][j] = a[i] * b[j];}}return U;}}void ArraSubs(vector< vector< double > > &C,vector< vector < double > > D){int a = C.size();int b = D.size();int c = C[0].size();int d = D[0].size();if (a!=b || c!=d){cerr << "Vectors not match in size !";}else{for (int i=1; i<a; i++){for (int j=1; j<a; j++){C[i][j] -= D[i][j];}}}}vector< double > VectSubs(vector< double > &a,vector< double > &b){int s1 = a.size();int s2 = b.size();if(s1 != s2){cerr << "Vectors not match in size !";}else{vector< double > value(s1,0);for (int i=1;i<s1;i++){value[i] = a[i] - b[i];}return value;}}void GausElim(vector< vector< double > > a) //高斯消元{int sz = a.size();vector< int > fx, ufx,ufxp, record;vector< double > x(sz, 0);vector< vector< double > > ret;//*****消元过程********for (int k = 1; k < sz-1; k++){double max = a[k][k];int p=0;for (int i=k+1; i<sz; i++){if (abs(a[i][k]) > abs(max)){p =i;max = a[i][k];}}//选出主元行if (abs(max) >= epsion){if (p != 0){for (int j=k; j<sz; j++){double temp = a[k][j];a[k][j] = a[p][j];a[p][j] = temp;}//交换主元行}for (int i=k+1; i<sz; i++){double tt = a[i][k] / a[k][k];for (int j=k+1; j<sz; j++){a[i][j] = a[i][j] - tt * a[k][j];}}}// end of if (abs(max) >= epsion)}// end of for (int k = 1; k < sz; k++)for (int i=1; i<sz; i++){int p=0;for (int j=1; j<=i;j++){if(abs(a[j][i]) >= 10*epsion)//不为零{p=j;}}record.push_back(p);}for (int i=1; i<sz; i++){int p=0;for (int j=sz-2; j>=0; j--){if (record[j] == i){p=j + 1;}}if (p != 0){ufxp.push_back(p);ufx.push_back(i);}}for (int i=1; i<sz; i++){int p = 0;for (int j=0; j < ufxp.size(); j++){if (ufxp[j] == i){p = 1;}}if (p == 0){fx.push_back(i);}}//end of for (int i=1; i<sz; i++)int c = fx.size();//***************************************************** for (int i=0; i<c; i++){for (int j=0; j<c; j++){if (i == j){x[fx.at(j)] = 1.0 + epsion;}else{x[fx.at(j)] = 0;}}int b = ufxp.size();for (int s=b-1; s>=0; s--){double temp=0;for(int t=ufxp.at(s)+1; t<sz; t++){temp += x[t] * a[ufx.at(s)][t];}x[ufxp.at(s)] = (0 - temp) / a[ufx.at(s)][ufxp.at(s)];}ret.push_back(x);}//end of for (int i=0; i<c; i++)//******************************************************* int sz1 = ret.size();int sz2 = ret[0].size();ofstream result2("CharactVector .txt", ios::app);result2.precision(12);for (int i=0; i<sz1; i++){for (int j=1;j<sz2; j++){result2 << scientific << ret[i][j] << endl;}result2 << endl << endl;}result2 << "下一个特征值的特征向量！" << endl;}void Stop(vector< vector< double > > &Ar){int a = Lambda.size();for (int i=0; i<a; i++){if (abs(Lambda[i].imag())<=epsion){vector< vector< double > > An=Ar;ArraSubs(An,NumbMultiArra(I,Lambda[i].real()));GausElim(An);}}ofstream result("result.txt");result.precision(12);vector<complex< double > >::iterator i,j;j = Lambda.end();for (i = Lambda.begin(); i<j; i++){result << scientific <<(*i) << endl;}exit(0);}void SolutS1S2(complex< double > &s1, complex< double > &s2,vector< vector< double > > &A){double s = A[m-1][m-1] + A[m][m];double t = A[m-1][m-1] * A[m][m] - A[m][m-1] * A[m-1][m];double det = s *s - 4.0 * t;if (det > 0){s1.imag(0);s1.real ((s + sqrt(det)) / 2);s2.imag (0);s2.real((s - sqrt(det)) / 2);}else{s1.imag(sqrt(0 - det) / 2);s1.real(s / 2);s2.imag(0 - sqrt(0 - det) / 2);s2.real(s / 2);}}void Save2(complex< double > &s1, complex< double > &s2)//存储特征值s1，s2 {Lambda.push_back(s1);Lambda.push_back(s2);}void Save1(complex< double > &s){Lambda.push_back(s);}void JudgemBelow2(vector< vector< double > > &A,vector< vector< double > > Abk){if (m == 1){complex< double > lbd;lbd.imag(0);lbd.real(A[1][1]);Save1(lbd);Stop(Abk);}else if (m == 0){Stop(Abk);}}//拟上三角化void Hessenberg(vector< vector< double > > &A){int a = A.size();vector< double > u(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;for (int r=1; r<a-2; r++){if (!IsBirZeroH(A, r)){d = GetModeH(A,r);c = (A[r+1][r] > 0) ?(0-d) : d;h = c * c - c * A[r+1][r];for (int j=1; j<a; j++){if (j < r+1){u[j] = 0;}else if (j == r+1){u[j] = A[j][r] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)p = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);q = NumbMultiVect(ArraMultiVect(A, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(A, VectMultiCovVect(w, u));ArraSubs(A, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("NISHANGDANJIAO.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}}void QRMethod(vector< vector< double > > A){int a = A.size();vector< vector< double > > Q(a,vector < double > (a,0));for (int i=1; i<a; i++){Q[i][i] = 1.0;}vector< vector< double > > RQ(A);//vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);vector< double > l(a,0);double d, c, h, t;for (int r=1; r<a-1; r++){if (!IsBirZero(A, r)){d = GetMode(A,r);c = (A[r][r] > 0) ?(0-d) : d;h = c * c - c * A[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = A[j][j] - c;}else{u[j] = A[j][r];}}// end of for (int j=1; j<=m; j++)l = ArraMultiVect(Q, u);ArraSubs(Q, VectMultiCovVect(l,NumbMultiVect(u, 1/h)));v = NumbMultiVect(ConveArraMultiVect(A, u), 1 / h);ArraSubs(A, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(RQ, u), 1 / h);q = NumbMultiVect(ArraMultiVect(RQ, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(RQ, VectMultiCovVect(w, u));ArraSubs(RQ, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)ofstream result("QR.txt");result.precision(12);for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << Q[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << A[i][j] << " ";}result << endl;}result << endl << endl;for (int i=1; i<a; i++){for (int j=1; j<a; j++){result << scientific << RQ[i][j] << " ";}result << endl;}}void CalculatAk(vector< vector < double > > &Ak){int a = Ak.size();vector< vector < double > > B(a,vector < double > (a,0));vector< double > u(a,0);vector< double > v(a,0);vector< double > p(a,0);vector< double > q(a,0);vector< double > w(a,0);double d, c, h, t;B = GetM(Ak);for (int r=1; r<a-1; r++){if (!IsBirZero(B, r)){d = GetMode(B,r);c = (B[r][r] > 0) ?(0-d) : d;h = c * c - c * B[r][r];for (int j=1; j<a; j++){if (j < r){u[j] = 0;}else if (j == r){u[j] = B[j][j] - c;}else{u[j] = B[j][r];}}// end of for (int j=1; j<=m; j++)v = NumbMultiVect(ConveArraMultiVect(B, u), 1 / h);ArraSubs(B, VectMultiCovVect(u, v));p = NumbMultiVect(ConveArraMultiVect(Ak, u), 1 / h);q = NumbMultiVect(ArraMultiVect(Ak, u), 1 / h);t = VectMultiVect(p, u) / h;w = VectSubs(q, NumbMultiVect(u, t));ArraSubs(Ak, VectMultiCovVect(w, u));ArraSubs(Ak, VectMultiCovVect(u, p));}// end of if (!IsBirZero(B, r))}// end of for (int r=1; r<m; r++)}int _tmain(int argc, _TCHAR* argv[]){vector< vector< double > > A(11, vector< double >(11,0));vector< vector< double > > Abk;I=A;for (int i=1; i<11; i++){vector< double > temp;for (int j=1; j<11; j++){if(i != j){A[i][j] = sin(0.5 * i + 0.2 * j);I[i][j] = 0;}else{A[i][j] = 1.5 * cos(i + 1.2 *j);I[i][j] = 1.0;}}}Abk = A;n = A.size() - 1;m = n;//初始化问题Hessenberg(A);QRMethod(A);while (1){if (abs(A[m][m-1]) <= epsion){complex< double > lbdm;lbdm.imag(0);lbdm.real(A[m][m]);Save1(lbdm);m -= 1;//对A进行降维处理！！！！A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();}JudgemBelow2(A, Abk);}else{complex< double > va1, va2;SolutS1S2(va1, va2, A);if (m == 2){Save2(va1,va2);Stop(Abk);}//end of if (m == 2)if ( abs(A[m-1][m-2]) <= epsion){Save2(va1,va2);m = m - 2;//矩阵降维A.pop_back();A.pop_back();int a = A.size();for (int i=0; i<a; i++){A[i].pop_back();A[i].pop_back();}JudgemBelow2(A, Abk);}else{if (k == L){cerr << "Stop without solution";exit(-1);}else{CalculatAk(A);k += 1;}}}// end of if (abs(A[m][m-1]) >= epsion)}}三、实验结果（1）A经过拟上三角化后得到的矩阵-8.827516758830e-001 -9.933136491826e-002 -1.103349285994e+000-7.600443585637e-001 1.549101079914e-001 -1.946591862872e+000-8.782436382928e-002 -9.255889387184e-001 6.032599440534e-0011.518860956469e-001-2.347878362416e+000 2.372370104937e+000 1.819290822208e+0003.237804101546e-001 2.205798440320e-001 2.102692662546e+0001.816138086098e-001 1.278839089990e+000 -6.380578124405e-001-4.154075603804e-0011.0556********e-016 1.728274599968e+000 -1.171467642785e+000-1.243839262699e+000 -6.399758341743e-001 -2.002833079037e+0002.924947206124e-001 -6.412830068395e-001 9.783997621285e-0022.557763574160e-001-5.393383812774e-017 -3.165873865181e-017 -1.291669534130e+000-1.111603513396e+000 1.171346824096e+000 -1.307356030021e+0001.803699177750e-001 -4.246385358369e-001 7.988955239304e-0021.608819928069e-0011.533464996622e-017 5.963321406181e-017 0.000000000000e+0001.560126298527e+000 8.125049397515e-001 4.421756832923e-001-3.588616128137e-002 4.691742313671e-001 -2.736595050092e-001 -7.359334657750e-0021.300562737229e-016 -3.097060010889e-017 0.000000000000e+0000.000000000000e+000 -7.707773755194e-001 -1.583051425742e+000-3.042843176799e-001 2.528712446035e-001 -6.709925401449e-0012.544619929082e-0011.610216724767e-016 -2.211571837369e-016 0.000000000000e+0000.000000000000e+000 6.483933100712e-017 -7.463453456938e-001-2.708365157019e-002 -9.486521893682e-001 1.195871081495e-0011.929265617952e-0021.368550186199e-016 7.151513190805e-017 0.000000000000e+0000.000000000000e+000 -2.522454384246e-017 1.072074718358e-016-7.701801374364e-001 -4.697623990618e-001 4.988259468008e-0011.137691603776e-001-2.780851300718e-017 -6.708630788363e-017 0.000000000000e+0000.000000000000e+000 -3.282796821065e-017 4.879774287475e-0171.854829286220e-016 7.013167092107e-001 1.582180688475e-0013.862594614233e-001-2.124604440055e-017 -1.707979758930e-016 0.000000000000e+0000.000000000000e+000 1.013496750890e-016 -4.153758325241e-0171.222621691887e-016 -5.551115123126e-017 4.843807602783e-0013.992777995177e-001（2）Q-3.519262579534e-001 4.427591982245e-001 -6.955982513607e-0016.486200753651e-002 3.709718861896e-001 1.855847143605e-001-1.628942319628e-002 -1.181053169648e-001 -5.255375383724e-002 -5.486582943568e-002-9.360277287361e-001 -1.664679186545e-001 2.615299548560e-001 -2.438671728934e-002 -1.394774360893e-001 -6.977585391241e-0026.124472142963e-003 4.440505443139e-002 1.975907909728e-0022.062836970533e-002-4.208697095111e-017 -8.810520554692e-001 -3.989762796959e-0013.720308728479e-002 2.127794064090e-001 1.064463557221e-001-9.343171079758e-003 -6.774200464527e-002 -3.014340698675e-002 -3.146955080444e-002-2.150178169911e-017 4.009681353156e-017 -5.371806806439e-001 -1.234945854205e-001 -7.063151608719e-001 -3.533456368496e-0013.101438948264e-002 2.248676491598e-001 1.000601783527e-0011.044622748702e-0016.113458775639e-018 -3.721194648197e-0177.068175586628e-0189.892235468621e-001 -1.239411731211e-001 -6.200358589825e-0025.442272839461e-003 3.945881637235e-002 1.755813350011e-0021.833059462907e-0025.184948268593e-017 -4.198303559531e-017 3.255071298412e-017-1.527665297025e-017 5.323610690264e-001 -6.733900344896e-0015.910581205868e-002 4.285425323867e-001 1.906901343193e-0011.990794495295e-0016.419444583601e-017 4.121668945107e-017 3.096964582901e-017-5.050360323651e-017 -7.078737686239e-017 -6.0597********e-001 -9.165783032818e-002 -6.645586508974e-001 -2.957110877580e-001 -3.0872********e-0015.455993559780e-017 -9.724896332186e-017 3.956603694870e-0171.913857001710e-018 -4.316583456405e-0172.797376665557e-0179.933396625117e-001 -9.690440311939e-002 -4.311990584470e-002-4.501694411183e-002-1.108640877071e-017 4.655237056115e-017 -1.0722********e-017 -9.470236121745e-018 4.277993227986e-017 8.866601870176e-017 -2.516725033065e-016 5.410088006061e-001 -5.817540838226e-001 -6.0734********e-001-8.470152033092e-018 9.650816729410e-017 -1.429362211392e-017 -2.789222052040e-017 -3.690793770141e-017 -8.090765462521e-017 -1.964050295697e-016 -6.772274683437e-017 -7.221591336735e-0016.917269588876e-001（3）R2.508342744917e+000 -2.185646885493e+000 -1.314609070786e+000-3.558787493836e-002 -2.609857850388e-001 -1.283121847090e+000 -1.390878610606e-001 -8.712897972161e-001 3.849367902971e-0013.353802899665e-0012.100627753398e-016 -1.961603277854e+000 2.407523727633e-0017.054714572823e-001 5.957204318279e-001 5.526978774676e-001-3.268209924413e-001 -5.769498668364e-002 2.871129330189e-001 -8.895128754189e-002-3.300197935770e-016 -1.872873410421e-016 2.404534601993e+0001.706758096328e+000 -4.239566704091e-001 3.405332305815e+000-1.050017655852e-001 1.462257102734e+000 -6.684487469283e-001 -4.027*********e-0013.0773********e-017 1.746386351950e-017 0.000000000000e+0001.577122080722e+000 6.399535133956e-001 3.468127872427e-001-5.701786649768e-002 4.014788054433e-001 -2.222476176311e-001 -6.317059236442e-0021.760039865880e-016 9.988285339980e-017 0.000000000000e+0000.000000000000e+000 -1.447846997770e+000 -1.415724007744e+000-2.806139044665e-001 -2.817910521892e-001 -4.611434881851e-0021.996629079956e-0018.804885435596e-017 4.996802050991e-017 0.000000000000e+0000.000000000000e+000 8.831633340975e-017 1.231641451542e+0001.619701003419e-001 1.962638275504e-001 5.350035621760e-001-1.509273424767e-001-7.728357669400e-018 -4.385868928757e-018 0.000000000000e+0000.000000000000e+000 -7.751835246841e-018 -1.279231078565e-017-7.753441914209e-001 -3.464514508821e-001 4.312226803504e-0011.234643696237e-001-5.603391361152e-017 -3.179943413314e-017 0.000000000000e+0000.000000000000e+000 -5.620413613517e-017 -9.274974944455e-0170.000000000000e+000 1.296312940612e+000 -4.288053318338e-0012.737334158165e-001-2.493361499851e-017 -1.414991023727e-017 0.000000000000e+0000.000000000000e+000 -2.500935953597e-017 -4.127119443934e-0170.000000000000e+000 0.000000000000e+000 -6.707396440648e-001-4.842320121884e-001-2.603055667460e-017 -1.477242832192e-017 0.000000000000e+0000.000000000000e+000 -2.610963355436e-017 -4.308689959101e-0170.000000000000e+000 0.000000000000e+000 -1.110223024625e-0167.168323926323e-002（4）RQ1.163074414164e+0002.632670934508e+000 -1.772796003272e+000-8.668899138521e-002 3.300503471047e-001 1.455162371214e+000 -9.730650448593e-001 -4.873031174655e-001 -7.756411630489e-001 -3.249201979113e-0011.836115060851e+000 1.144286420080e-001 -9.880381403133e-0015.589725694767e-001 4.694190067101e-002 -2.978478237007e-0011.617130577649e-002 6.936977702522e-001 1.367670571405e-0011.419099231519e-0025.598200524418e-016 -2.118520153533e+000 -1.876189745783e+000-5.407071940597e-001 1.171538359721e+000 -2.550323020223e+0001.691577936540e+000 1.229951613262e+000 1.387947777212e+0008.667502917242e-001-3.179059303049e-017 -5.261714527977e-017 -8.471995127808e-0014.382910468318e-001 -1.008632199185e+000 -7.959374261495e-0014.769258865577e-001 4.072683083890e-001 4.096390493527e-0013.363378940862e-001-3.514195999269e-016 3.391949386044e-017 -2.160938214545e-016 -1.432244342447e+000 -5.742284908055e-001 1.213151477723e+000 -3.457508625575e-001 -4.749853573124e-001 -3.176158274191e-001 -4.294507015032e-002-3.704735750656e-017 -1.0998********e-016 -1.953241363386e-0178.269089741494e-017 6.556779598004e-001 -9.275250974463e-0012.529079844053e-001 6.905949216976e-001 -2.374430675823e-002-2.429781119781e-001-6.420051822287e-017 3.865817713597e-017 -3.806718328740e-0172.129734126613e-017 7.853*********e-017 4.698400884876e-001-2.730776009527e-001 7.821296259798e-001 -9.580964936399e-0027.846239841323e-0021.478496020372e-016 -8.383952267131e-017 9.577540382744e-017-7.120338911078e-017 -1.220876056602e-016 -1.854471832043e-0161.287679058937e+000 -3.576058900348e-001 -4.116725408807e-0033.914268216423e-0014.477158378610e-017 -6.204017427568e-017 3.360322998507e-017-1.133882337819e-017 -2.759056827456e-017 -9.217582575819e-0172.214639994444e-016 -3.628760503545e-001 7.398980975354e-0017.241608309576e-0023.408891482791e-017 2.353495494255e-017 1.932283926688e-017-3.456910153081e-017 -2.015177337156e-017 -8.084422691634e-017 -5.839476980893e-017 5.551115123126e-017 -5.176670596524e-0024.958522909877e-002（5）特征值(6.360627875745e-001,0.000000000000e+000)(5.650488993501e-002,0.000000000000e+000)(9.355889078188e-001,0.000000000000e+000)(-9.805309562902e-001,1.139489127430e-001)(-9.805309562902e-001,-1.139489127430e-001)(1.577548557113e+000,0.000000000000e+000)(-1.484039822259e+000,0.000000000000e+000)(-2.323496210212e+000,8.930405177200e-001)(-2.323496210212e+000,-8.930405177200e-001)(3.383039617436e+000,0.000000000000e+000)（6）实特征值的特征向量4.745031936539e+0003.157868541750e+0001.729946912420e+001-1.980049331455e+000-3.187521973524e+0017.794009603201e+000-1.004255685845e+0011.670757770480e+0011.310524273054e+0011.000000000001e+000下一个实特征值对应的特征向量：-5.105003830625e+000-4.886319842360e+0009.505161576151e+000-6.788331788214e-001-9.604334110499e+000-3.0457********e+0001.574873885602e+001-7.395037126277e+000-7.109654943661e+0001.000000000001e+000下一个实特征值对应的特征向量：2.792418944529e+0001.598236841512e+000-5.207507040911e-001-1.667886451702e+000-1.225708535859e+0017.241214790777e+000-5.398214501433e+0002.841008912974e+001-1.216518754416e+0011.000000000001e+000下一个实特征值对应的特征向量：6.217350824581e-001-1.115111815226e-001-2.483773580804e+000-1.306860840421e+000-3.815605442533e+0008.117305509395e+000-1.239170883675e+000-6.800309586197e-0012.691900144840e+0001.000000000001e+000下一个实特征值对应的特征向量：-1.730784582112e+0012.408192534967e+0014.133852365119e-001-8.572044074538e+0009.287334657928e-002-7.832726042776e-002-6.374274025716e-001-3.403204760832e-001。
西南交大数值分析第二次大作业(可以运行)
数值分析第二次大作业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、矩阵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;}三、计算结果。
数值分析第二次作业答案
练习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。
数值分析-第二次作业答案
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 。
数值分析第二次大作业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。
北航数值分析报告大作业二
数值分析大作业（二）学院名称宇航学院专业名称航空宇航推进理论与工程学生姓名段毓学号SY16153062016年11月5日1 算法设计方案首先将矩阵A 进行拟上三角化，把矩阵A 进行QR 分解，计算出RQ 。
要得出矩阵A 的全部特征值，首先对A 进行QR 的双步位移得出特征值。
最后，采用列主元的高斯消元法求解特征向量。
1.1 A 的拟上三角化因为对矩阵进行QR 分解并不改变矩阵的结构，因此在进行QR 分解前对矩阵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 Λ执行 若()n r r i a r ir,,3,2)(Λ++=全为零，则令)()1(r r A A =+,转5；否则转2。
计算()∑+==nri r ir r 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 +-=令()nTr nrr r r r r r r r R a a c a u ∈-=++)()(,2)(,1,,,,0,,0ΛΛ。
计算r r T r r h u A p /)(=r r rr 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(继续。
1.2 A 的QR 分解具体算法如下所示，记)1(1-=n A A ，并记[]nn r ij r a A ⨯=)(，令I Q =1 对于1,,2,1-=n r Λ执行 1.若()n r r i a r ir ,,3,1)(Λ++=全为零，则令r r Q Q =+1r r A A =+1,转5；否则转2。
北航 数值分析第二次大作业(带双步位移的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 分解，构造如下形式的矩阵：
清华大学高等数值分析 第二次实验作业
400
600
800
1000
1200
1400
1600
1800
2000
迭代次数
图4
m步的重启动GMRES法求解Ax=b的收敛曲线
结论： m步重启GMRES方法快于m步重启Arnoldi方法， 随m增加， 迭代次数减小， 但都大于未重启算法的次数。当m=40时两种方法计算时间最短，此外，m步重 启动 Arnoldi 方法的收敛曲线有峰点和波纹，收敛速度均匀性差， m 步重启动 GMRES方法的收敛曲线很平滑，单调下降，收敛速度均匀性好。（图4是五条曲 线， 只是由于m=20和m=80两条曲线比较靠近， 看起来像四条， 放大后才能看清） 2.对于 1 中的矩阵，将特征值进行平移，使得实部有正有负，和 1 的结果进行比 较，方法的收敛速度会如何？基本的 Arnoldi 算法有无峰点？若有，基本的 GMRES 算法相应地会怎样？ 解： （1）欲将特征值进行平移，使得实部有正有负，可以将矩阵 A 做如下变换：
10
0
特征值虚部按不同比例因子 k变化的 GMRES算法的收敛曲线 (阶数 n=1000）
k=0.2 k=0.5 k=2 k=5
10
-5
||rk||/||b||
10
-10
10
-15
0
100
200
300
400
500
600
700
800
900
1000
迭代次数
图8 特征值虚部按不同比例因子k变化的GMRES法求解求解Ax=b的收敛曲线
图7 特征值虚部按不同比例因子k变化的Arnoldi法求解求解Ax=b的收敛曲线
（3）GMRES法求解求解A′x=b: 利用matlab编程实现GMRES算法。b = ones(1000,1)，x0 = zeros(1000, 1)。 计算每一步迭代的残差rk相对于初始残差的2范数。相对残差2范数的对数 值与迭代步数的关系曲线如图8所示：
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
《数值分析》计算实习报告第二题院系：机械工程及自动化学院学号：姓名：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 ++=-∈。
（4）计算()()(1)()///r T r r r r r r rTr r r rr r r rr r T T r r r rp A u h q A u h t p u h q t u A A u u p ωω+====-=--（5）继续。
当此算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n−1)。
2.用带双步位移的QR 方法求矩阵A 全部特征值为了加速收敛，采用带双步位移的QR 方法求解A 的全部特征值，具体算法如下：（1）由步骤1得到拟上三角矩阵(1)n A -；给定精度水平0ε> 和最大迭代次数L 。
（2）记(1)(1)1=n ij n nA A a -⨯⎡⎤=⎣⎦，令k = 1,m = n 。
（3）如果(),1k m m a ε-≤，则得到A 的一个特征值()k mm a ，置m:=m-1（降阶），转（4）；否则转（5）。
（4）如果2m ≤，则得到矩阵A 的一个特征值()11k a ，或两个特征值12,s s ，转（9）；否则转（3）。
（5）如果()1,2k m m a ε--≤，则得到A 的两个特征值12,s s ，置m : = m – 2（降阶），转（4）；否则转（6）。
（6）如果k = L ，则终止计算，未得到A 的全部特征值；否则转（7）。
（7）记()(1,)k k ij m m A a i j m ⨯⎡⎤=≤≤⎣⎦，计算()()1,1()()()()1,1,11,21)() (k k m m mmk k k k m m mm m m m mk k k k k k k T k k k ks a a t a a a a M A sA tI I m M Q R M QR A Q A Q ------+=+=-=-+==是阶单位矩阵对作分解 （8）置k : = k + 1，转（3）。
（9）A 的全部特征值已计算完毕，停止计算。
其中，k M 的QR 分解与1k A +的计算用下列算法实现：记(1)()11[],[],r k ij m m r ij m m k B M b B b C A ⨯⨯====。
对于1,2,,1r m =-执行（1）若()()1,2,,r ir b i r r m =++全为零，则令11,r r r r B B C C ++==，转（5）；否则转（2）。
（2）计算()()()()1,2()sgn 0,r r r r rr r r r r r r r r r rrd c b d a c d h c c b +==-===-若则取 （3）令()()()()1,0,,0,,,,Tr r r m r rr r r r mr u b c b b R +=-∈。
（4）计算11////T r r r r T r r r r T r r r r r r r rTr r r rr r r rT T r r r r r rv B u h B B u v p C u h q C u h t p u h q t u C C u u p ωω++==-====-=--（5）继续。
此算法执行完后，得到1k m A C +=。
3.用列主元素Gauss 消去法求特征向量记λ为矩阵A 的实特征值，x 为对应的特征向量，则A x =λx ，即(A - λ)x = 0的解即为矩阵A 的特征向量。
因此对于特征向量的求解可采用列主元素Gauss 消去法。
0A I λ-=，经检查由于矩阵A 无重特征值，所以每个特征值对应的特征向量的线性空间为1，矩阵A I λ-的秩为N-1，经过列主元素Gauss 消去法消元后，矩阵有且仅有最后一行全为0。
此时设定1n x =-，再依次求出()1,2,,1k x k n n =--，就可得到对应λ的其中一个特征向量12(,,,)T n x x x x =，对应λ的全部特征向量可表示为k x 。
最后为了得到标准形式，将向量x 正交化。
具体算法如下。
1. 消元过程 对于1,2,,1k n =-执行（1）选行号k i ，使()()max k k k i k ik k i na a ≤≤=。
（2）交换()k kj a 与()(,1,,)k k i j a j k k n =+所含的数值。
（3）对于1,2,,i k k n =++计算()()(1)()()/ (1,2,,)k k ik ik kkk k k ijij ik kim a a aam a j k k n +==-=++2. 回代过程()()()111,2,,1n n k k k kj j kk j k x x a x a k n n =+=-⎛⎫=-=-- ⎪⎝⎭∑最终得到的向量12(,,,1)T n x x x =-即为A 对应于实特征值λ的特征向量。
3. 将向量正交化。
三、源程序代码#include<windows.h>#include<iostream> #include<iomanip> #include<math.h> using namespace std;#define N 10 #define E 1e-12#define MAX 10000 //最大迭代次数//创建原始矩阵Avoid create_A(double a[N][N]) { int i, j; for (i = 0; i<N; i++) for (j = 0; j<N; j++) if (i != j) a[i][j] = sin(0.5*(i + 1) + 0.2*(j + 1)); else a[i][j] = 1.52*cos((i + 1) + 1.2*(j + 1)); }//符号函数 int sgn(double a) { int z; if (a>E) z = 1; else12(,,,1)T n x x x =-z = -1;return z;}//矩阵中小于或等于E的元素置0void zhi0(double A[N][N]){int i, j;for (i = 0; i<N; i++)for (j = 0; j<N; j++)if (fabs(A[i][j]) <= E)A[i][j] = 0;}//二次多项式求deltadouble delta(double b, double c){double m;m = b*b - 4 * c;return m;}//求二次方程两个根s1，s2void s1s2(double A[N][N], int m, double lamdaR[N], double lamdaI[N]) {double s, t, x;s = 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];x = delta(s, t);if (x>E){x = sqrt(x);lamdaR[m] = s / 2 - x / 2;lamdaR[m - 1] = s / 2 + x / 2;lamdaI[m] = 0;lamdaI[m - 1] = 0;}else{x = sqrt(fabs(x));lamdaR[m] = s / 2;lamdaR[m - 1] = s / 2;lamdaI[m] = -x / 2;lamdaI[m - 1] = x / 2;}}//矩阵的拟上三角化void nishangsjh(double A[N][N]){int i, j, k, flag;double d, c, h, t;double u[N], p[N], q[N], w[N];for (i = 0; i<N - 2; i++){flag = 0;for (j = i + 2; j<N; j++)if (fabs(A[j][i])>E){flag = 1;break;}if (flag == 0)continue;for (j = i + 1, d = 0; j<N; j++)d = d + A[j][i] * A[j][i];d = sqrt(d);c = -1 * sgn(A[i + 1][i])*d;h = c*c - c*A[i + 1][i];for (j = i + 2; j<N; j++)u[j] = A[j][i];for (j = 0; j<i + 2; j++)u[j] = 0;u[i + 1] = A[i + 1][i] - c;for (j = 0; j<N; j++){for (k = i + 1, p[j] = 0; k<N; k++)p[j] = A[k][j] * u[k] + p[j];p[j] = p[j] / h;}for (j = 0; j<N; j++){for (k = i + 1, q[j] = 0; k<N; k++)q[j] = A[j][k] * u[k] + q[j];q[j] = q[j] / h;}for (j = 0, t = 0; j<N; j++)t = t + p[j] * u[j];t = t / h;w[j] = q[j] - t*u[j];for (j = 0; j<N; j++)for (k = 0; k<N; k++)A[j][k] = A[j][k] - w[j] * u[k] - u[j] * p[k];}zhi0(A);}//对矩阵A进行一次基本QR方法迭代void QRdiv(double A[N][N],double RQ[N][N]){int i, j, k, flag;double Q[N][N], R[N][N];double d, c, h;double u[N], w[N], p[N];for (i = 0; i<N; i++)for (j = 0; j<N; j++)if (i == j)Q[i][j] = 1;elseQ[i][j] = 0;for (i = 0; i<N; i++)for (j = 0; j<N; j++)R[i][j] = A[i][j];for (i = 0; i<N - 1; i++){flag = 0;for (j = i + 1; j<N; j++)if (fabs(R[j][i])>E){flag = 1;break;}if (flag == 0)continue;for (j = i, d = 0; j<N; j++)d = d + R[j][i] * R[j][i];d = sqrt(d);c = -1 * sgn(R[i][i])*d;h = c*c - c*R[i][i];for (j = i + 1; j<N; j++)u[j] = R[j][i];u[j] = 0;u[i] = R[i][i] - c;for (j = 0; j<N; j++)for (k = 0, w[j] = 0; k<N; k++)w[j] = Q[j][k] * u[k] + w[j];for (j = 0; j<N; j++)for (k = 0; k<N; k++)Q[j][k] = Q[j][k] - w[j] * u[k] / h;for (j = 0; j<N; j++){for (k = i, p[j] = 0; k<N; k++)p[j] = R[k][j] * u[k] + p[j];p[j] = p[j] / h;}for (j = 0; j<N; j++)for (k = 0; k<N; k++)R[j][k] = R[j][k] - u[j] * p[k];}zhi0(R);zhi0(Q);for (i = 0; i<N; i++)for (j = 0; j<N; j++)for (k = 0, RQ[i][j] = 0; k<N; k++)RQ[i][j] = RQ[i][j] + R[i][k] * Q[k][j]; }//矩阵M的QR分解void QRdiv_M(double M[N][N], double A[N][N], int &m) {int i, j, k, flag;double d, c, h, tr;double u[N], w[N], p[N], q[N], v[N];for (i = 0; i<m; i++){flag = 0;for (j = i + 1; j<N; j++)if (fabs(A[j][i])>E){flag = 1;break;}if (flag == 0)continue;for (j = i, d = 0; j <= m; j++)d = d + M[j][i] * M[j][i];d = sqrt(d);c = -1 * sgn(M[i][i])*d;h = c*c - c*M[i][i];for (j = i + 1; j <= m; j++)u[j] = M[j][i];for (j = 0; j<i; j++)u[j] = 0;u[i] = M[i][i] - c;for (j = 0; j <= m; j++){for (v[j] = 0, k = i; k <= m; k++)v[j] = v[j] + M[k][j] * u[k];v[j] = v[j] / h;}for (k = i; k <= m; k++)for (j = 0; j <= m; j++)M[k][j] = M[k][j] - u[k] * v[j];for (j = 0; j <= m; j++){for (p[j] = 0, k = i; k <= m; k++)p[j] = p[j] + A[k][j] * u[k];p[j] = p[j] / h;}for (j = 0; j <= m; j++){for (q[j] = 0, k = i; k <= m; k++)q[j] = q[j] + A[j][k] * u[k];q[j] = q[j] / h;}for (tr = 0, j = i; j <= m; j++)tr = tr + p[j] * u[j];tr = tr / h;for (j = 0; j <= m; j++)w[j] = q[j] - tr*u[j];for (j = 0; j <= m; j++)for (k = 0; k <= m; k++)A[j][k] = A[j][k] - (w[j] * u[k] + u[j] * p[k]);}}//带双步位移的QR方法求特征值int tezhengzhi(double A[N][N], double lamdaR[N], double lamdaI[N]) {int L = 0, m = N - 1;int i, j, k;double s, t;double M[N][N];int flag;loop3:if (fabs(A[m][m - 1]) <= E){lamdaR[m] = A[m][m]; lamdaI[m] = 0;m--;goto loop4;}elsegoto loop5;loop4:if (m<2){if (m == 0){lamdaR[0] = A[0][0];lamdaI[0] = 0;}if (m == 1)s1s2(A, m, lamdaR, lamdaI);flag = 1;goto loop9;}elsegoto loop3;loop5:if (fabs(A[m - 1][m - 2]) <= E){s1s2(A, m, lamdaR, lamdaI);m = m - 2;goto loop4;}elsegoto loop6;loop6:if (L == MAX){flag = 0;goto loop9;}elsegoto loop7;loop7:s = 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];for (i = 0; i <= m; i++)for (j = 0; j <= m; j++)if (i == j){for (k = 0, M[i][j] = 0; k <= m; k++)M[i][j] = A[i][k] * A[k][j] + M[i][j];M[i][j] = M[i][j] - s*A[i][j] + t;}else{for (k = 0, M[i][j] = 0; k <= m; k++)M[i][j] = A[i][k] * A[k][j] + M[i][j];M[i][j] = M[i][j] - s*A[i][j];}QRdiv_M(M, A, m);goto loop8;loop8:L++;goto loop3;loop9:return L;}//列主元的高斯消元法求解特征向量void gauss(double x[N][N], double tzR[N], int &nR, double lamdaR[N], double lamdaI[N]) {double A[N][N];double sum, ar, ch, m[N];int i, j, k, p, ik;for (i = 0, nR = 0; i<N; i++)if (fabs(lamdaI[i]) <= E)tzR[nR++] = lamdaR[i];for (p = 0; p<nR; p++){create_A(A);for (i = 0; i<N; i++)for (j = 0; j<N; j++)if (i == j)A[i][j] = A[i][j] - tzR[p];for (k = 0; k<N - 1; k++){for (ik = k, i = k + 1; i<N; i++)if (A[ik][k]<A[i][k])ik = i;for (j = k; j<N; j++){ch = A[ik][j];A[ik][j] = A[k][j];A[k][j] = ch;}for (i = k + 1; i<N; i++){m[i] = A[i][k] / A[k][k];for (j = k + 1; j<N; j++)A[i][j] = A[i][j] - m[i] * A[k][j];}}x[p][N - 1] = -1;for (k = N - 2; k >= 0; k--){for (ar = 0, j = k + 1; j<N; j++)ar = ar + A[k][j] * x[p][j];x[p][k] = -ar / A[k][k];}for (j = 0, sum = 0; j<N; j++)sum = sum + x[p][j] * x[p][j];sum = sqrt(sum);for (j = 0; j<N; j++)x[p][j] = x[p][j] / sum;}}//打印N*N阶矩阵void print_A(double a[][N]){int i, j;printf("\n");for (i = 0; i<N; i++){for (j = 0; j<int(N / 2); j++)cout << setw(22) << a[i][j];}printf("\n");for (i = 0; i<N; i++){for (j = int(N / 2); j<N; j++)cout << setw(22) << a[i][j];printf("\n");}printf("\n");}//打印矩阵所有特征值void print_tzz(double lamdaR[N], double lamdaI[N]){int i;printf("\n");for (i = 0; i<N; i++){cout << "\tλ" << setw(2) << i + 1 << " = ";cout << setw(20) << lamdaR[i];if (fabs(lamdaI[i])>E)if (lamdaI[i]>E)cout << " + " << setw(20) << lamdaI[i] << setw(2) << "i";elsecout << " - " << setw(20) << -lamdaI[i] << setw(2) << "i";printf("\n");}printf("\n");}//打印实特征值的特征向量void print_tzxl(double x[N][N], double tzR[N], int nR){int i, j;for (i = 0; i<nR; i++){cout << endl << " 实特征值：" << setw(25) << tzR[i] << endl;cout << " 特征向量：" << endl;for (j = 0; j<int(N / 2); j++)cout << setw(22) << x[i][j];printf("\n");for (j = int(N / 2); j<N; j++)cout << setw(22) << x[i][j];}}int main(){double A[N][N], RQ[N][N], x[N][N];double lamdaR[N], lamdaI[N], tzR[N];int flag, nR;system("Color f0");cout << setiosflags(ios::scientific) << setiosflags(ios::right) << setprecision(12) << endl;create_A(A);cout << " 矩阵A的拟上三角化矩阵A(n-1)：" << endl;nishangsjh(A);print_A(A);printf("\n");cout << " 对矩阵A(n-1)进行QR分解，矩阵R*Q：" << endl;QRdiv(A, RQ);print_A(RQ);printf("\n");cout << " 矩阵A(n-1)经过QR迭代得到的最终矩阵：" << endl;flag = tezhengzhi(A, lamdaR, lamdaI);print_A(A);printf("\n");if (flag == 0)cout << "未求出所有特征值!" << endl;else{cout << " 矩阵A的全部特征值：" << endl;print_tzz(lamdaR, lamdaI);printf("\n");cout << " 实特征值对应的特征向量：" << endl;gauss(x, tzR, nR, lamdaR, lamdaI);print_tzxl(x, tzR, nR);}return 0;}四、计算结果1. 拟上三角矩阵A(n−1)如下：2. 对矩阵A(n−1)进行一次基本QR方法迭代，得到如下矩阵：3. 矩阵A(n−1)经过QR迭代得到的最终矩阵如下：4. 矩阵A的全部特征值如下：5.对应于实特征值的特征向量如下：五、结果分析与讨论1. 由计算结果1和结果2可知，对拟上三角矩阵A(n−1)进行一次基本QR方法迭代，所得矩阵仍然为拟上三角阵，这与两个矩阵相似的理论推导结果相符，证明了拟上三角化算法的正确性。