李津高等数值分析计算实践第一次大作业

合集下载

数值分析大作业一

数值分析大作业一一、算法设计方案1、求λ1和λ501的值：思路：采用幂法求出按模最大特征值λmax，该值必为λ1或λ501，若λmax小于0，则λmax=λ1；否则λmax=λ501。

再经过原点平移，使用幂法迭代出矩阵A-λmax I的特征值，此时求出的按模最大特征值即为λ1和λ501的另一个值。

2、求λs的值：采用反幂法求出按模最小的特征值λmin即为λs，其中的方程组采用LU分解法进行求解。

3、求与μk最接近的特征值：对矩阵A采用带原点平移的反幂法求解最小特征值，其中平移量为：μk。

4、A的条件数cond(A)=| λmax/λmin|；5、A的行列式的值：先将A进行LU分解，再求U矩阵对角元素的乘积即为A 行列式的值。

二、源程序#include<iostream>#include<iomanip>#include<math.h>#define N 501#define E 1.0e-12 //定义精度常量#define r 2#define s 2using namespace std;double a[N];double cc[5][N];void init();double mifa();double fmifa();int max(int aa,int bb);int min(int aa,int bb);int max_3(int aa,int bb,int cc);void LU();void main(){double a1,a2,d1,d501=0,ds,det=1,miu[39],lamta,cond;int i,k;init();/*************求λ1和λ501********************/a1=mifa();if(a1<0)d1=a1; //若小于0则表示λ1的值elsed501=a1; //若大于0则表示λ501的值for(i=0;i<N;i++)a[i]=a[i]-a1;a2=mifa()+a1;if(a2<0)d1=a2; //若小于0则表示λ1的值elsed501=a2; //若大于0则表示λ501的值cout<<"λ1="<<setiosflags(ios::scientific)<<setprecision(12)<<d1<<"\t";cout<<"λ501="<<setiosflags(ios::scientific)<<setprecision(12)<<d501<<endl;/**************求λs*****************/init();ds=fmifa();cout<<"λs="<<setiosflags(ios::scientific)<<setprecision(12)<<ds<<endl;/**************求与μk最接近的特征值λik**************/cout<<"与μk最接近的特征值λik:"<<endl;for(k=0;k<39;k++){miu[k]=d1+(k+1)*(d501-d1)/40;init();for(i=0;i<N;i++)a[i]=a[i]-miu[k];lamta=fmifa()+miu[k];cout<<"λi"<<k+1<<"\t\t"<<setiosflags(ios::scientific)<<setprecision(12)<<lamta<<en dl;}/**************求A的条件数**************/cout<<"矩阵A的条件式";cond=abs(max(abs(d1),abs(d501))/ds);cout<<"cond="<<setiosflags(ios::scientific)<<setprecision(12)<<cond<<endl;/**************求A的行列式**************/cout<<"矩阵A的行列式";init();LU();for(i=0;i<N;i++){det*=cc[2][i];}cout<<"det="<<setiosflags(ios::scientific)<<setprecision(12)<<det<<endl;system("pause");}/**************初始化函数，给a[N]赋值*************/void init(){int i;for(i=1;i<=501;i++)a[i-1]=(1.64-0.024*i)*sin((double)(0.2*i))-0.64*exp((double)(0.1/i)); }/**************幂法求最大绝对特征值**************/double mifa(){int i,k=0;double u[N],y[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++) //控制最大迭代次数为2000{/***求y(k-1)***/double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;}/****求新的uk****/u[0]=a[0]*y[0]+b*y[1]+c*y[2];u[1]=b*y[0]+a[1]*y[1]+b*y[2]+c*y[3]; //前两列和最后两列单独拿出来求中D间的循环求for(i=2;i<N-2;i++){u[i]=c*y[i-2]+b*y[i-1]+a[i]*y[i]+b*y[i+1]+c*y[i+2];}u[N-2]=c*y[N-4]+b*y[N-3]+a[N-2]*y[N-2]+b*y[N-1];u[N-1]=c*y[N-3]+b*y[N-2]+a[N-1]*y[N-1];/***求beta***/double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}//cout<<"Beta"<<k<<"="<<Beta<<"\t"; 输出每次迭代的beta /***求误差***/error=abs(Beta-Beta_)/abs(Beta);if(error<=E) //若迭代误差在精度水平内则可以停止迭代{return Beta;} //控制显示位数Beta_=Beta; //第个eta的值都要保存下来，为了与后个值进行误差计算 }if(k==2000){cout<<"error"<<endl;return 0;} //若在最大迭代次数范围内都不能满足精度要求说明不收敛}/**************反幂法求最小绝对特¬征值**************/double fmifa(){int i,k,t;double u[N],y[N]={0},yy[N]={0},b=0.16,c=-0.064,Beta_=0,error;for(i=0;i<501;i++)u[i]=1; //令u[N]=1for(k=1;k<2000;k++){double sum_u=0,gh_sum_u;for(i=0;i<N;i++){sum_u+=u[i]*u[i]; }gh_sum_u=sqrt(sum_u);for(i=0;i<N;i++){y[i]=u[i]/gh_sum_u;yy[i]=y[i]; //用重新赋值，避免求解方程组的时候改变y的值}/****LU分解法解方程组Au=y，求新的***/LU();for(i=2;i<=N;i++){double temp_b=0;for(t=max(1,i-r);t<=i-1;t++)temp_b+=cc[i-t+s][t-1]*yy[t-1];yy[i-1]=yy[i-1]-temp_b;}u[N-1]=yy[N-1]/cc[s][N-1];for(i=N-1;i>=1;i--){double temp_u=0;for(t=i+1;t<=min(i+s,N);t++)temp_u+=cc[i-t+s][t-1]*u[t-1];u[i-1]=(yy[i-1]-temp_u)/cc[s][i-1];}double Beta=0;for(i=0;i<N;i++){Beta+=y[i]*u[i];}error=abs(Beta-Beta_)/abs(Beta);if(error<=E){return (1/Beta);}Beta_=Beta;}if(k==2000){cout<<"error"<<endl;return 0;} }/**************求两数最大值的子程序**************/int max(int aa,int bb){return(aa>bb?aa:bb);}/**************求两数最小值的子程序**************/int min(int aa,int bb){return(aa<bb?aa:bb);}/**************求三数最大值的子程序**************/int max_3(int aa,int bb,int cc){ int tt;if(aa>bb)tt=aa;else tt=bb;if(tt<cc) tt=cc;return(tt);}/**************LU分解**************/void LU(){int i,j,k,t;double b=0.16,c=-0.064;/**赋值压缩后矩阵cc[5][501]**/for(i=2;i<N;i++)cc[0][i]=c;for(i=1;i<N;i++)cc[1][i]=b;for(i=0;i<N;i++)cc[2][i]=a[i];for(i=0;i<N-1;i++)cc[3][i]=b;for(i=0;i<N-2;i++)cc[4][i]=c;for(k=1;k<=N;k++){for(j=k;j<=min(k+s,N);j++){double temp=0;for(t=max_3(1,k-r,j-s);t<=k-1;t++)temp+=cc[k-t+s][t-1]*cc[t-j+s][j-1];cc[k-j+s][j-1]=cc[k-j+s][j-1]-temp;}//if(k<500){for(i=k+1;i<=min(k+r,N);i++){double temp2=0;for(t=max_3(1,i-r,k-s);t<=k-1;t++)temp2+=cc[i-t+s][t-1]*cc[t-k+s][k-1];cc[i-k+s][k-1]=(cc[i-k+s][k-1]-temp2)/cc[s][k-1];}}}}三、程序结果。

JZX高等数值分析第一次实验作业

相对残差 6.1302e-16 8.7797e-09 8.0295e-09 8.5677e-09 9.1433e-09
a、 m=1 (左为相对残差，右为取对数情况)
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
0 1
1.5
0
-5
-10
-15
-20
-25
-30
-35
-40
1000
4.761361
2.9675e-11
5.7069e-04
3.5336e+11
1500
15.646571
1.0778e-11
6.8236e-04
2.3672e+10
3000
132.198513
6.3164e-13
1.1700e-04
2.8110e+11
a、n=1000 时步数与相对残差关系图：（上为相对残差，下为取对数情况下结果）
（2）当 A 最大特征值远大于第二个特征值，最小特征值远小于第二个最小特征值时收敛
性情况。
思路：构造题目要求的矩阵 A。首先随机生成 n 阶矩阵 B，B 不满秩，构造对角阵 A1（最
大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值），则由此构
造出对称正定矩阵 A： b1=B’*B; A=b1’*A1*b1。同样设定精确解 Xj 为元素全部为 1 的 n
5、构造对称不定的矩阵，验证 Lanczos 方法的近似中断，观察收敛曲线中的峰点个数和特
征值的分布关系；观察当出现峰点时，MINRES 方法的收敛性态怎样。
解：思路：类似前两题，首先构造出一个 n 阶对角阵 D，其对角线上有 m 个负值，再对随

清华大学高等数值计算(李津)实践题目二(SVD计算及图像压缩)(包含matlab代码)

方法3：方法
主要思想：
传统的算法和分而治之方法都是基于先将矩阵二对角化的方法，而方法是完全不同的方法，它通过一系列平面旋转（也称为变换）；最终使得矩阵收敛于一个对角矩阵。
其中每一个设计成，使得对于一对选定的指标，有如下式子成立：
最终收敛时，令，有，且具有互相正交的列向量，可以将写成，其中为对角矩阵，为正交矩阵；从而得到原矩阵的分解式。
先将矩阵二对角化得到二对角阵，然后将求的问题分为两个子问题。
先将矩阵写成如下形式：
其中，为上二对角矩阵，是相应维数的向量中的第个单位向量，并且一般我们取。
现在假如我们已经求得了，的如下：
，
并且令为的最后一行，为的第一行。那么将这些带入，可以得到：
可以看出中间的矩阵形式上很简单，除了对角元素和第行上的元素，其它元素都为；这里先把它的的计算问题放到最后，而假定已知其为：，将它带入上式，就可以得到的分解式：
主要步骤：
令是最初的矩阵，是经过一次旋转过得到的矩阵，其中为对角矩阵，是列向量范数全为1的矩阵。可以知道对进行单边的变换相当于对矩阵进行双边的变换。最终得到的矩阵的列向量的范数就是所要求的奇异值，并且最终归一化以后的矩阵的列向量就是左奇异向量。
第2部分实验计算
针对以上三种算法用进行数值模拟，三个算法源自工具箱，分别为（传统的SVD算法），（分而治之）和（方法）。实验中找10幅彩色图片（尺寸为），并将其灰度化，转化为数据矩阵，计算3种算法的计算时间和最小奇异值误差，以及左右奇异向量误差，并得到了相应的对比曲线图。左右奇异向量的误差通过和的矩阵2-范数的来评估。最小特征值误差，使用默认自带的SVD算法得到的结果作为参考标准（实际中大多数情况也是无法知道确切的值）。

数值分析第一次作业

《数值分析》计算作业院系：航空科学与工程学院学号： SY1005512姓名：王天龙日期： 2010年10月31日计算实习说明书目的：训练运用计算机进行科学与工程计算的能力。

要求：1．独立进行算法设计、程序设计和上机运算，并得出正确的结果。

2．编制程序时全部采用双精度，要求按题目的要求设计输出，并执行打印。

3．只能根据题目给出的信息并且只允许一次计算得出全部结果。

题目：第一题设有501×501的矩阵123499500501a b c b a b cc b a b c A c b a b c c b a b c ba ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦其中0.1(1.640.024)s i n (0.2)0.64 (125i i a i i e i =--= ，，，；0.16b =；0.064c =-。

矩阵A 的特征值12501()ii λ= ，，，满足 125011501||min ||S i i λλλλλ≤≤<<<= ，试求：1.1λ，501λ和S λ的值。

2.A 的与数5011140k kλλμλ-=+最接近的特征值(1239)ik k λ= ，，，。

3.A 的（谱范数）条件数2()cond A 和行列式det A 。

说明：1．在所有的算法中，凡是要给出精度水平ε的，都取1210ε-=。

2．选择算法时，应使A 的所有零元素都不存储。

3．打印以下内容：（1）算法的设计方案。

（2）全部源程序（要求注明主程序和每个子程序的功能）。

（3）特征值1λ，501λ，S λ和(1239)ik k λ= ，，，以及2()cond A ，det A 的值。

4.采用e 型输出所有计算结果，并至少显示12位有效数字。

一、程序算法的设计算法设计方案如下：二、全部源程序编程软件：Fortran：三、计算结果1.特征值1λ，501λ，S λ1-.107001135582E+02λ=，501 .972463398616E+01λ=，-.555823879237E-02S λ=2. (1239)ik k λ= ，，，如下表所示（ZK 代表ik λ）3.A 的条件数2()cond A 和行列式det A 的值2() .192509100000E+04cond A =，det .277059968428+119A =五、讨论这里选取的初始向量为X(i)=1,X={x1,x2,x3,,,,,,,x501},当初始向量与特征向量较近时，收敛较快，若初始向量与特征向量正交，则求解可能失真。

数值分析计算实习作业一

数值分析计算实习题一学号：：院系：2015年11月5日一、分析1.1算法分析题目要求求出：1）特征值从小到大排列的最小特征值1λ和最大特征值501λ。

2）特征值中模最小的特征值s λ。

3）靠近一组数k μ的一组特征值k i λ。

4）矩阵A 的条件数cond(A)2。

5）行列式detA 。

解决方法：1）若将所有行列式按模的大小排列则模最大的特征值一定是1λ和501λ中的一个，因此利用幂法求出模最大的特征值1m λ。

然后利用带原点平移的幂法，将系数矩阵变为1m A I λ-即将所有特征值都减去1m λ，则特征值按大小顺序排列的次序不变，模最大的特征值依然在整个排列的两端，再用一次幂法得到模最大的特征值21=m m λλλ-，其中λ为带原点平移的幂法求出的特征值，最后两个特征值1m λ、2m λ比较大小，大的为501λ，小的为1λ。

2）因为s λ为按模最小的特征值，因此用反幂法可求的其特征值。

3）因为k i λ靠近数k μ，因此k i k λμ-一定是所有的k λμ-中模最小的，因此可利用带原点平移的反幂法求出特征值k i λ，此时的系数矩阵变为k A I μ-。

4）条件数cond(A)2为模最小的特征值与模最大的特征值的比的绝对值，因此利用1和2中求出的1m λ和s λ可解出条件数。

5）可对矩阵A 进行LU 分解，即A LU =则det()det()det()A L U =⨯，又因为矩阵L 对角线元素为1，则det()L =1，所以det()det()A U =，U 为上三角阵，行列式为对角线元素的乘积，因此可得A 的行列式。

1.2程序分析1.2.1 因为A 为拟三角阵，储存时零元素不储存，因此将矩阵A 压缩为5*501的矩阵CA 的带元素ij a =C 中的元素1,i j s j c -++ 程序中A[5][501]即为压缩后的矩阵。

1.1.2 程序中的B[5][501]为过渡矩阵，在幂法迭代、反幂法迭代以及LU 分解中均用矩阵B 来计算，计算之间对B 进行适当的赋值。

清华大学高等数值分析作业李津1——矩阵基础

20130917题目求证：在矩阵的LU 分解中，111n n Tn ij i j j i j L I e e α-==+⎛⎫=- ⎪⎝⎭∑∑证明：在高斯消去过程中，假设0jj a ≠ ，若a=0，可以通过列变换使得前面的条件成立，这里不考虑这种情况。

对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=••-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。

()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。

故111n n T n ij i j n j i j L I e e I α-==+⎛⎫⎛⎫=- ⎪ ⎪ ⎪⎝⎭⎝⎭∏∑上式中的每一项均是初等变换，从右向左乘，则每乘一次相当于对右边的矩阵进行一次向下乘法叠加的初等变换。

由于最初的矩阵为单位阵，变换从右向左展开，因而每一次变换不改变已经更新的数据，既该变换是从右向左一列一列更新数据，故11nn Tn ij i j j i j L I e e α==+⎛⎫=- ⎪⎝⎭∑∑。

数学证明：1n Tij i j i j e e α=+⎛⎫ ⎪⎝⎭∑具有,000n j j A -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j A B --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此： 1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……，，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。

(完整版)数值分析第一次作业

问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868；（2）S ’’(0.25)=S ’’(0.53)=0。

分析：本问题是已知五个点，由这五个点求一三次样条插值函数。

边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。

对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。

⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]）解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。

数值分析第一次大作业

《数值分析》计算实习报告第一题院系：机械工程及自动化学院_学号： _____姓名： _ ______2017年11月7日一、算法设计方案1、求λ1，λ501和λs 的值1)利用幂法计算出矩阵A 按模最大的特征值，设其为λm 。

2)令矩阵B =A −λm I （I 为单位矩阵），同样利用幂法计算出矩阵B 按模最大的特征值λm ′。

3)令λm ′′=λm ′+λm 。

由计算过程可知λm 和λm ′′分别为矩阵A 所有特征值按大小排序后，序列两端的值。

即，λ1=min⁡{λm ,λm ′′}，λ501=max⁡{λm ,λm ′′}。

4) 利用反幂法计算λs 。

其中，反幂法每迭代一次都要求解线性方程组1k k Au y -=，由于矩阵A 为带状矩阵，故可用三角分解法解带状线性方程组的方法求解得到k u 。

2、求A 的与数μk =λ1+k λ501−λ140最接近的特征值λi k (k =1,2, (39)1) 令矩阵D k =A −μk I ，利用反幂法计算出矩阵D k 按模最小的特征值λi k ′,则λi k =λi k ′+μk 。

3、求A 的（谱范数）条件数cond(A )2和行列式det A1) cond(A)2=|λm λs |，前文已算出m λ和s λ，直接带入即可。

2) 反幂法计算λs 时，已经对矩阵A 进行过Doolittle 分解，得到A=LU 。

而L 为对角线上元素全为1的下三角矩阵，U 为上三角矩阵，可知det 1L =，5011det ii i U u ==∏，即有5011det det det ii i A L U u ====∏。

最后，为节省存储量，需对矩阵A 进行压缩，将A 中带内元素存储为数组C [5][501]。

二、源程序代码#include<windows.h>#include<iostream>#include<iomanip>#include<math.h>using namespace std;#define N 501#define K 39#define r 2#define s 2#define EPSI 1.0e-12//求两个整数中的最大值int int_max2(int a, int b){return(a>b ? a : b);}//求两个整数中的最小值int int_min2(int a, int b){return(a<b ? a : b);}//求三个整数中的最大值int int_max3(int a, int b, int c){int t;if (a>b)t = a;else t = b;if (t<c) t = c;return(t);}//定义向量内积double dianji(double x[], double y[]) {double sum = 0;for (int i = 0; i<N; i++)sum = sum + x[i] * y[i];return(sum);}//计算两个数之间的相对误差double erro(double lamd0, double lamd1){double e, d, l;e = fabs(lamd1 - lamd0);d = fabs(lamd1);l = e / d;return(l);}//矩阵A的压缩存储初始化成Cvoid init_c(double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)if (i == 0 || i == 4)c[i][j] = -0.064;else if (i == 1 || i == 3)c[i][j] = 0.16;elsec[i][j] = (1.64 - 0.024*(j + 1))*sin(0.2*(j + 1)) - 0.64*exp(0.1 / (j + 1)); }//矩阵复制void fuzhi_c(double c_const[][N], double c[][N]){int i, j;for (i = 0; i<r + s + 1; i++)for (j = 0; j<N; j++)c[i][j] = c_const[i][j];}//LU三角分解void LUDet_c(double c_const[][N], double c_LU[][N]){double sum;int k, i, j;fuzhi_c(c_const, c_LU);for (k = 1; k <= N; k++){for (j = k; j <= int_min2(k + s, N); j++){sum = 0;for (i = int_max3(1, k - r, j - s); i <= k - 1; i++)sum += c_LU[k - i + s][i - 1] * c_LU[i - j + s][j - 1];c_LU[k - j + s][j - 1] -= sum;}for (j = k + 1; j <= int_min2(k + r, N); j++){sum = 0;for (i = int_max3(1, j - r, k - s); i <= k - 1; i++)sum += c_LU[j - i + s][i - 1] * c_LU[i - k + s][k - 1];c_LU[j - k + s][k - 1] = (c_LU[j - k + s][k - 1] - sum) / c_LU[s][k - 1];}}}//三角分解法解带状线性方程组void jiefc(double c_const[][N], double b_const[], double x[]){int i, j;double b[N], c_LU[r + s + 1][N], sum;for (i = 0; i<N; i++)b[i] = b_const[i];LUDet_c(c_const, c_LU);for (i = 2; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= i - 1; j++)sum += c_LU[i - j + 2][j - 1] * b[j - 1];b[i - 1] -= sum;}x[N - 1] = b[N - 1] / c_LU[2][N - 1];for (i = N - 1; i >= 1; i--){sum = 0;for (j = i + 1; j <= int_min2(i + 2, N); j++)sum += c_LU[i - j + 2][j - 1] * x[j - 1];x[i - 1] = (b[i - 1] - sum) / c_LU[2][i - 1];}}//幂法求按模最大特征值double mifa_c(double c_const[][N]){double u[N], y[N];double sum, length_u, beta0, beta1;int i, j;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;for (i = 1; i <= N; i++){sum = 0;for (j = int_max2(i - 2, 1); j <= int_min2(i + 2, N); j++)sum = sum + c_const[i - j + 2][j - 1] * y[j - 1];u[i - 1] = sum;}beta1 = dianji(u, y);} while (erro(beta0, beta1) >= EPSI);return(beta1);}//反幂法求按模最小特征值double fmifa_c(double c_const[][N]){double u[N], y[N];double length_u, beta0, beta1;int i;for (i = 0; i<N; i++)//迭代初始向量u[i] = 0.5;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);do{beta0 = beta1;length_u = sqrt(dianji(u, u));for (i = 0; i<N; i++)y[i] = u[i] / length_u;jiefc(c_const, y, u);beta1 = dianji(y, u);} while (erro(beta0, beta1) >= EPSI);beta1 = 1 / beta1;return(beta1);}//计算lamd_1、lamd_501、lamd_svoid calculate1(double c_const[][N], double &lamd_1, double &lamd_501, double &lamd_s) {int i;double lamd_mifa0, lamd_mifa1, c[r + s + 1][N];lamd_mifa0 = mifa_c(c_const);fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - lamd_mifa0;lamd_mifa1 = mifa_c(c) + lamd_mifa0;if (lamd_mifa0<lamd_mifa1){lamd_1 = lamd_mifa0;lamd_501 = lamd_mifa1;}else{lamd_501 = lamd_mifa0;lamd_1 = lamd_mifa1;}lamd_s = fmifa_c(c_const);}//平移+反幂法求最接近u_k的特征值void calculate2(double c_const[][N], double lamd_1, double lamd_501, double lamd_k[]){int i, k;double c[r + s + 1][N], h, temp;temp = (lamd_501 - lamd_1) / 40;for (k = 1; k <= K; k++){h = lamd_1 + k*temp;fuzhi_c(c_const, c);for (i = 0; i<N; i++)c[2][i] = c[2][i] - h;lamd_k[k - 1] = fmifa_c(c) + h;}}//计算cond(A)和det(A)void calculate3(double c_const[][N], double lamd_1, double lamd_501, double lamd_s, double &cond_A, double &det_A){int i;double c_LU[r + s + 1][N];if (fabs(lamd_1)>fabs(lamd_501))cond_A = fabs(lamd_1 / lamd_s);elsecond_A = fabs(lamd_501 / lamd_s);LUDet_c(c_const, c_LU);det_A = 1;for (i = 0; i<N; i++)det_A *= c_LU[2][i];}//*主程序*//int main(){int i, count = 0;double c_const[5][N], lamd_k[K];double lamd_1, lamd_501, lamd_s;double cond_A, det_A;//设置白背景黑字system("Color f0");//矩阵A压缩存储到c[5][501]init_c(c_const);cout << setiosflags(ios::scientific) << setiosflags(ios::right) << setprecision(12) << endl;//计算lamd_1、lamd_501、lamd_scalculate1(c_const, lamd_1, lamd_501, lamd_s);cout << " 矩阵A的最小特征值：λ1 = " << setw(20) << lamd_1 << endl;cout << " 矩阵A的最大特征值：λ501 = " << setw(20) << lamd_501 << endl;cout << " 矩阵A的按模最小的特征值：λs = " << setw(20) << lamd_s << endl;//求最接近u_k的特征值calculate2(c_const, lamd_1, lamd_501, lamd_k);cout << endl << " 与数u_k最接近的特征值：" << endl;for (i = 0; i<K; i++){cout << " λ_ik_" << setw(2) << i + 1 << " = " << setw(20) << lamd_k[i] << " ";count++;if (count == 2){cout << endl;count = 0;}}//计算cond_A和det_Acalculate3(c_const, lamd_1, lamd_501, lamd_s, cond_A, det_A);cout << endl << endl;cout << " 矩阵A的条件数：cond(A) = " << setw(20) << cond_A << endl;cout << " 矩阵A的行列式的值：det(A) = " << setw(20) << det_A << endl << endl;return 0;}三，计算结果四，分析初始向量选择对计算结果的影响当选取初始向量0(1,1,,1)Tu=时，计算的结果如下：此结果即为上文中的正确计算结果。

李津高等数值分析计算实践第一次大作业

一、生成矩阵1.5个100阶对角方阵。

根据下述问题，设计5个100阶对角矩阵，对角元分布如下：表1.1 5个对角阵的对角元分布生成的5个对角矩阵中，前四个矩阵为正定矩阵，第五个矩阵为非正定矩阵。

其中，前三个矩阵条件数都为100。

2.生成10个100阶矩阵生成10个100阶矩阵，做QR分解，生成10个Q矩阵。

由于生成的Q,即对应的特征向量，对最后的求解结果影响不大，故此处随机生成了10个100阶矩阵，且为了每次的结果一样，将矩阵存到文件中。

一、实验结果（1）用共轭梯度法，Lanczos法，MINRES法进行计算，分析特征值分布对算法的影响。

计算得到的结果如下表所示：其中的为e w最后一次迭代之后求得的解与真实值的差，即e w=x w−x exact.从表格中可以看出以下几点：●在矩阵特征值相同时，特征值的分布对收敛性有影响。

在特征值偏向较小特征值的时候，CG法和Minres算法迭代次数增加，而Lanczos算法更是出现了不收敛的情况。

相对的，当特征值偏向较大特征值的时候，三种算法的迭代次数都均匀分布的有所下降。

故可有结论，特征值偏向最大值时有较好收敛性。

●矩阵条件数对收敛性有影响。

由三种算法的第3，4组数据能看出，3对应矩阵的条件数大于4矩阵的条件数，且3的迭代次数均小于4的迭代次数，故可有结论，条件数大的正定矩阵有较好的收敛性。

●矩阵是否正定对收敛性有影响。

由于CG法不适用于非正定矩阵，故只对Minres法和Lanczos法分析。

Lanczos法在实验中的非正定矩阵的运算中不收敛，且Minres法的迭代次数也较之前有了大幅度的提高。

故可有结论，正定矩阵较非正定矩阵有较好的收敛性。

（2）取定一个对角矩阵，用10个单位对称阵做实验，观察特征向量对算法的影响。

此处为了观察特征向量的影响，选定了最普通的1号对角阵。

结果如下表所示：表2.2 共轭梯度法对于不同特征向量的结果表2.4 Minres法对于不同特征向量的结果从以上三个表格数据可以看出，矩阵对于不同特征向量的变化收敛性变化不大，在10个不同特征向量下，迭代次数变化最大在3次之内，可以认为对迭代次数没有影响。

清华大学高等数值计算(李津)实践题目一(共轭梯度CG法_Lanczos算法与MINRES算法)

高等数值计算实践题目一1. 实践目的本次计算实践主要是在掌握共轭梯度法，Lanczos 算法与MINRES 算法的基础上，进一步探讨这3种算法的数值性质，主要研究特征值特征向量对算法收敛性的影响。

2. 实践过程（一）生成矩阵（1）作5个100阶对角阵i D 如下：1D 对角元：1,1,...,20,1+0.1(-20),21,...,100j j d j d j j ====2D 对角元：1,1,...,20,1+(-20),21,...,100j j d j d j j ==== 3D 对角元：,1,...,80,81,81,...,100j j d j j d j ====4D 对角元：,1,...,40,41,41,...,60,41+(60),61,...,100j j j d j j d j d j j =====-= 5D 对角元：,1,...,100j d j j ==记i D 的最大模特征值和最小模特征值分别为1iλ和in λ，则i D 特征值分布有如下特点：1D 的特征值有较多接近于i n λ，并且1/i i n λλ较小，2D 的特征值有较多接近于i n λ，并且1/i i n λλ较大， 3D 的特征值有较多接近于1i λ，并且1/i i n λλ较大，4D 的特征值有较多接近于中间模特征值，并且1/i i n λλ较大， 5D 的特征值均匀分布，并且1/i i n λλ较大（2）随机生成10个100阶矩阵j M ：(100(100))j M fix rand =g并作它们的QR 分解，得j Q 和j R ，这样可得50个对称的矩阵Tij j i j A Q DQ =，其中i D 的对角元就是ij A 的特征值，若它们都大于0，则ij A 正定，j Q 的列就是相应的特征向量。

结合（1）可知，ij A 都是对称正定阵。

（二）计算结果以下计算，均选定精确解(100,1)exact x ones =，初值0(100,1)x zeros =由ij exact kA x b =计算得到k b （算法中要求解的精度为10e -）。

数值分析计算实习第一题

直接用定义: ��(��)2 = ‖��‖2‖��−1‖2
求 A 的条件数很繁琐,需要先进行化简：
首先：
由于 A 是对称矩阵，
‖��‖2 = ��max(��)
说明：
1. 在所用的算法中，凡是要给出精度水平的ε，都取 ��=10−12。
2. 选择算法的时候应使矩阵 A 的所有零元素都不存储。
3. 打印以下内容：
(1)算法设计方案和思路。
(2)全部源程序。
(3)特征值��1,��501,��,��(��=1,2,⋯,39)以及��(��)2, det��的值(采用 e 型输出实型数，并至少显示 12 位有效数字)。
λi[16] -2.533970311130E+00 λi[38] 8.648666065193E+00
λi[17] -2.003230769563E+00 λi[39] 9.254200344575E+00
λi[18] -1.503557611227E+00 cond(A)2 1.925204273903E+03
λi[19] -9.935586060080E-01 det(A) 2.772786141752E+118

高等数值分析作业-第一次实验

高等数值分析作业-第一次实验-CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN-140-120-100-80-60-40-20020迭代次数l g （|r k |）高等数值分析第一次实验T1. 构造例子说明CG 的数值形态。

当步数 = 阶数时CG 的解如何当A 的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何 Answer:对于问题1：当步数 = 阶数时CG 的解如何在MATLAB 中构造N 阶对称正定矩阵代码如下：N=1000D = diag(rand(N,1)); U = orth(rand(N,N)); A = U'*D*U;在计算时，取X0=zeros(N,1);b=ones(N,1);自己编写CG 算法，如下： Xk = X0; rk=b-A*Xk; pk=rk; crk_1=rk'*rk; for k=1:N k=k+1; apk=A*pk;ak=crk_1/(pk'*apk); Xk=Xk+ak*pk; rk=rk-ak*apk; crk=rk'*rk; bk_1=crk/crk_1; crk_1=crk; pk=rk+bk_1*pk; m(k)=norm(rk); r(k)=k; endplot(r,m,'r-'); Ek=m(k)计算结果如下（绘制出来的log 10‖‖随迭代次数的变化如上图所示）：N 1000 2000 3000 4000 log 10‖rr ‖ 运行时间（s ）由上表可以看出对于对称正定矩阵A ，CG 算法还是比较稳定的，但求解步数=阶数时，CG 算法的解即为准确解（误差极小）。

对于问题2：当A 的最大特征值远大于第二个最大特征值，最小特征值远小于第二个最小特征值时，方法的收敛性如何构造1000阶的对称正定矩阵如下，收敛准则取为绝对ε<10(-10)：首先构造一个特征值分别为到1的对称正定矩阵A ，代码如下（算例1）：D = diag(linspace,1,N)); U = orth(rand(N,N)); A = U'*D*U;在之前的基础上，将最大特征值调为107，最小特征值为10-5，代码如下（算例2）：DIA=linspace,1,N); DIA(1)=10^(-5); DIA(N)=10^7; D = diag(DIA); U = orth(rand(N,N)); A = U'*D*U;最后生成一个特征值在10-5到107均匀分布的矩阵（算例3）:DIA=linspace(10^(-5),10^7,N); D = diag(DIA); U = orth(rand(N,N)); A = U'*D*U;计算结果如右图所示，首先对比可以发现矩阵的收敛速度跟其条件数大小有关，条件数小时，收敛速度快，算例1>2>3，同时，A 的中间特征值分布对CG 的收敛速度有巨大的影响。

清华大学高等数值分析(李津)所有作业答案合集

对矩阵A 进行LU 分解，()()()()()1111111L M n M M M n ---=-=∙∙-………… ，其中()1n Tn ij i j i j M j I e e α=+⎛⎫=+ ⎪⎝⎭∑ ，i e 、j e 为n 维线性空间的自然基。

()M j 是通过对单位阵进行初等变换得到，通过逆向的变换则可以得到单位阵，由此很容易得到()M j 的逆矩阵为1n Tn ij i j i j I e e α=+⎛⎫- ⎪⎝⎭∑。

数学证明：1nTi j i ji j ee α=+⎛⎫ ⎪⎝⎭∑具有,000n j jA -⎛⎫ ⎪⎝⎭ 和1,1000n j n j B -+-+⎛⎫⎪⎝⎭ 的形式，且有+1,-11,10000=000n j j n j n j AB --+-+⎛⎫⎛⎫⎪⎪⎝⎭⎝⎭ 而11n n T ij i j j k i j e e α-==+⎛⎫ ⎪⎝⎭∑∑具有1,1000n k n k B -+-+⎛⎫⎪⎝⎭的形式，因此：1311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+⎡⎤⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫=---⎢⎥ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪⎝⎭⎝⎭⎢⎥⎝⎭⎝⎭⎝⎭⎝⎭⎣⎦⎛⎫⎛⎫⎛⎫=-- ⎪ ⎪ ⎝⎭⎝⎝⎭∏∑∏∑∑∑∑∑……11211n n n T Tk n ik i kk k i k e I e e α--===+⎛⎫⎛⎫=- ⎪⎪ ⎪⎭⎝⎭⎝⎭∑∑∑#20130924题目一问：能否用逐次householder 相似变换变实矩阵A 为上三角矩阵，为什么？解：不能用逐次householder 相似变换变A 为上三角矩阵，原因如下：A 记作：()12=,,n A a a a ……，，存在householder 阵1H s.t. 1111H a e α= ，则()()()111111111111111111111,,,0T Th H AH H a A H e H A H e H A H h H A H ααα⎛⎫'''=== ⎪⎪'⎝⎭⎛⎫''=+ ⎪ ⎪⎝⎭11H A H ''第一列的元素不能保证为1e 的倍数，故无法通过householder 变换实现上三角化。

数值分析第一次作业

数值分析第一次作业班级学号姓名习题24、用Newton法求方程f(x)=x^3-2*x^2-4*x-7=0在[3，4]中的根。

代码：function[x_star,k]=Newton1[fname,dfname,x0,ep,Nmax]if nargin<5 Nmax=500; endif nargin<4 ep=1e-5;endx=x0;x0=x+2*ep;k=0;while abs(x0-x)>ep&k<Nmax k=k+1x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);endx_star=x;if k==Nmax warning(‘已迭代上限次数’);endfname=inline('x^3-2*x^2-4*x-7');dfname=inline('3*x^2-4*x-4');[x_star,k]=Newton1(fname,dfname,3.5)x_star =3.6320k =4方法二：2－4用割线法求方程的根function [x_star,k]=Gline(fun,x0,x1,ep,Nmax)if nargin<5 Nmax=500;endif nargin<4 ep=1e-5;endk=0;while abs(x1-x0)>ep&k<Nmaxk=k+1;x2=x1-feval(fun,x1)*(x1-x0)/(feval(fun,x1)-feval(fun,x0))x0=x1;x1=x2;endx_star=x1;if k==Nmax warning('已迭代上限次数');endfun=inline('x^3-2*x^2-4*x-7');[x_star,k]=Gline(fun,3,4)x2 =3.5263x2 =3.6168x2 =3.6327x2 =3.6320x2 =3.6320x_star =3.6320k =5习题33、用列主元消去法解方程组[-1 2 -2; 3 -1 4; 2 -3 -2][x1 x2 x3]=[-1 7 0] 代码：function x=Gauss_x1(A,b)A=[A’;b]’,n=length(b);for k=1:n-1s=A(k,k);p=k;for i=l+1:nif abs(s)<abs(A(i,k))s=A(I,k);p=I;endendAfor i=k+1:nm=A(i,k)/A(k,k);fprintf(‘m%d%d=%f\n’,i,k,m);for j=k:n+1A(i,j)=A(i,j)-m*A(k,j);endendfprintf(‘A%d=\n’,k+1);AendA(n,n+1)=A(n,n+1)/A(n.n);for i=n-1:-1:1s=0for j=i+1:ns=s+A(i,j)*A(j,n+1);endA(i,n+1)=(A(i,n+1)-s)/A(i,i);endA(:,n+1)A=[-1,2,-2;3,-1,4;2,-3,-2];b=[-1;7;0];x=Gauss_x1(A,b)A =3.0000 -1.00004.0000 7.00000 1.6667 -0.6667 1.33330 -2.3333 -4.6667 -4.6667A=3.0000 -1.00004.0000 7.00000 -2.3333 -4.6667 -4.66670 0 -4.0000 -2.0000x =2.00001.00000.50004、用追赶法解三对角方程[2 -1 0 0 0;-1 2 -1 0 0;0 -1 2 -1 0;0 0 -1 2 -1;0 0 0 -1 2][x1 x2 x3 x4 x5=[1 0 0 0 0 ] 代码：function x=zhuigan(A,B,C,D)n=length(B);Xzeros(1,n);U=zeros(1,n);Q=zeros(1,n);U(1)=C(1)/B(1);Q(1)=D(1)/B(1);for i=2:n-1U(i)=C(i)/(B(i)-U(i-1)*A(i-1));endfor i=2:nQ(i)=(D(i)-Q(i-1)*A(i-1))/(B(i)-U(i-1)*A(i-1));endX(n)=Q(n);for i=n-1:-1:1X(i)=Q(i)-U(i)*X(i+1);endXA=[-1，-1，-1，-1;];B=[2,2,2,2,2];C=[-1,-1,-1,-1];D=[1;0;0;0;0];X=zhuigan(A,B,C,D)X= 0.8333 0.6667 0.5000 0.3333 0.16676、用三角分解法解方程组[-2 4 8;-4 18 -16;-6 2 -20][x1 x2 x3]=[5 8 7]代码function[y,x]=LU_s(A,b)b=b';A=[A';b]',n=length(b');x=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);endendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=b(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*x(j);endx(i)=(y(i)-ly1)/U(i,i);endA=[-2,4,8;-4,18,-16;-6,2,-20];b=[5;8;7];[y,x]=LU_s(A,b)A =-2 4 8 5-4 18 -16 8-6 2 -20 7L =1 0 02 1 03 -1 1U =-2 4 80 10 -320 0 -76y =5-2-10x =-1.53160.22110.13163-8用LU分解法解线性方程组[5,7,9,10;6,8,10,9;7,10,8,7;5,7,6,5][x1 x2 x3 x4]=[1 1 1 1] 代码function[y,x]=LU_s(A,b)b=b';A=[A';b]',n=length(b');x=zeros(n,1);y=zeros(n,1);U=zeros(n);L=eye(n);for k=1:nU(1,k)=A(1,k);L(k,1)=A(k,1)/U(1,1);endfor i=2:nfor k=i:nlu=0;lu1=0;for j=1:i-1lu=lu+L(i,j)*U(j,k);lu1=lu1+L(k,j)*U(j,i);endU(i,k)=A(i,k)-lu;L(k,i)=(A(k,i)-lu1)/U(i,i);EndendLUfor i=1:nly=0;for j=1:ily=ly+L(i,j)*y(j);endy(i)=b(i)-ly;endfor i=n:-1:1ly1=0;for j=i+1:nly1=ly1+U(i,j)*x(j);endx(i)=(y(i)-ly1)/U(i,i);endA=[5,7,9,10;6,8,10,9;7,10,8,7;5,7,6,5];b=[1;1;1;1];[y,x]=LU_s(A,b)A =5 7 9 10 16 8 10 9 17 10 8 7 15 76 5 1L =1.0000 0 0 01.2000 1.0000 0 01.4000 -0.5000 1.0000 01.0000 0 0.6000 1.0000U =5.0000 7.0000 9.0000 10.0000 0 -0.4000 -0.8000 -3.0000 0 0 -5.0000 -8.5000 0 0 0 0.1000y =1.0000 -0.2000 -0.5000 0.3000x =20.0000 -12.0000 -5.0000 3.0000。

829607103_2_高值第一次实验作业

高等数值分析实验作业一2021.11.02
(注意：矩阵规模不得小于1000 阶)
1.构造例子说明CG 的数值性态. 当步数= 阶数时CG 的解如何? 当A 的最
大特征值远大于第二个最大特征值, 最小特征值远小于第二个最小特征值时方法的收敛性如何?
2.对于同样的例子, 比较CG 和Lanczos 的计算结果.
3.当A 只有m 个不同特征值时, 对于大的m 和小的m, 观察有限精度
下Lanczos 方法如何收敛.
4.取初始近似解为零向量, 右端项b 仅由A 的m 个不同个特征向量的线性组
合表示时, Lanczos 方法的收敛性如何? 数值计算中方法的收敛性和m 的大小关系如何?
5.构造对称不定的矩阵, 验证Lanczos 方法的近似中断, 观察收敛曲线中的峰
点个数和特征值的分布关系; 观察当出现峰点时, MINRES 方法的收敛性态怎样.
1。

数值分析第一次作业解答

数值分析第一次作业解答1：(a) —个问题的病态性如何，与求解它的算法有关系。

x ;(b) 无论问题是否病态，好的算法都会得到它好的近似解。

x ；(C)计算中使用更高的精度，可以改善问题的病态性。

X ；(d) 用一个稳定的算法计算一个良态问题，一定会得到他好的近似解。

V;(e) 浮点数在整个数轴上是均匀分布。

x ;(f) 浮点数的加法满足结合律。

x(g) 浮点数的加法满足交换律。

X ;(h) 浮点数构成有效集合。

V;(i) 用一个收敛的算法计算一个良态问题，一定得到它好的近似解2: 程序t=0.1;n=1:10;e=n/10-n*te = 1.0e-015 *[ 0 0 -0.0555 0 0-0.1110 -0.1110 0 0 0] 由舍人误差造成n=3，6，7 时的结果不为零。

4:两种等价的一元二次方程求解公式-b - Pb2 - 4acx =2a2cx 二-b b2 - 4ac对a=1, b=-100000000, c=1，应采用哪种算法?A二[1,-100000000,1];roots(A);可得：X1 = 100000000;x2=0a=1;b=-100000000;c=1;x1仁(-b-sqrt(b*b-4*a*c))/(2*a)x12=(-b+sqrt(b*b-4*a*c))/(2*a)x2仁2*c/(-b-sqrt(b*b-4*a*c))x22=2*c/(-b+sqrt(b*b-4*a*c))由第一种算法：X1 = 100000000;x2=7.45058 X10由第二种算法：X1 = 13417728;x2=-1.0 X108原因：太小的数作分母。

5:程序:fun cti on y=tt(x)s=0;t=x;n=1;while s+t~=s;s=s+t;t=-x A2/(( n+1)*( n+2))*tn=n+2;endntt(2n 1)eps)(a)t小于计算机的计算精度。

李津 高等数值分析 计算实践第一次大作业

数值分析大作业一

JZX高等数值分析第一次实验作业

清华大学高等数值计算(李津)实践题目二(SVD计算及图像压缩)(包含matlab代码)

数值分析第一次作业

数值分析计算实习作业一

清华大学高等数值分析作业李津1——矩阵基础

(完整版)数值分析第一次作业

数值分析第一次大作业

李津高等数值分析计算实践第一次大作业

清华大学高等数值计算(李津)实践题目一(共轭梯度CG法_Lanczos算法与MINRES算法)

数值分析计算实习第一题

高等数值分析作业-第一次实验

清华大学高等数值分析(李津)所有作业答案合集

数值分析第一次作业

829607103_2_高值第一次实验作业

数值分析第一次作业解答

李津高等数值分析计算实践第一次大作业