《数值分析B》
计算实习题目第二题
院系：013交通科学与工程学院姓名：康乐
学号：SY0913107
一、算法设计方案
1、按照题目要求将矩阵初始化；
2、按照带双步位移的QR 方法求实矩阵特征值的算法步骤，计算目标矩阵的特征值。 （1）将目标矩阵拟上三角化。编写拟上三角化子程序，在主函数中直接调用。拟上三角化子程序编写算法：
①对于一个10*10的a 矩阵，先设定循环次数为8次；
②第一次循环选中a 矩阵第一列，并将第一列第一行元素取0，其他元素不变，即
102090[0,,,.......]T a a a ，然后计算所有元素平方和s ，取c ，根据10a 的符号判断若10a 是正值
则c= -s ，若10a 是负值或者0，则c=s ，然后设定u[10],将u[1]元素减去c 值，其他的u[i]（i=0，2……9）元素不变；设定H[10][10],H=2/()T T I uu u u -；进而计算(2)
a
HaH =，即a 矩阵
左右各乘以H 矩阵然后再赋给a 矩阵。
③第二次是选择2步求出来的a 矩阵的第二列，然后依次按照2的步骤循环八次，最终得到的a 矩阵就是拟上三角矩阵。
（2）输出拟上三角化以后的矩阵a ；
（3）进入带双步位移的QR 方法循环，其中要用到QR 分解方法，故先给出QR 分解子程序，再给出具体的带双步位移的QR 方法子程序。
QR 分解子程序：
①对于10*10的a 矩阵，先设定循环次数为九次，将需QR 分解的矩阵HK[10][10]初始化为单位矩阵，该矩阵保存最后的变换矩阵；
②第一次循环选中a 矩阵第一列，即00102090[,,,.......]T a a a a ，然后计算所有元素平方和s ，取c ，根据10a 的符号判断若10a 是正值则c= -s ，若10a 是负值或者0，则c=s ，然后设定u[10],将u[0]元素减去c 值，其他的u[i]（i=1，2……9）元素不变；设定H[10][10],H=2/()T
T
I uu u u -；进而计算(2)
a
HaH =，即a 矩阵左右各乘以H 矩阵然后再
赋给a 矩阵。
③第二次是选择2步求出来的a 矩阵的第二列，然后依次按照2的步骤循环九次，最终得到的a 矩阵就是上三角矩阵。
带双步位移的QR 方法子程序：
①先设置循环的初始条件m=9，m>=0；
②然后分情况讨论：当m>1的时候，又分为三种情况，第一种是a[m][m-1]的绝对值小于E 的时候，可以得出第m 个特征值为a[m][m],第二种情况是
|a[m][m-1]|>E,|a[m-1][m-2]| )() [ ]m m m m m m m m a a a a --- -组成的矩阵的两个特征值即为第 m ，m-1的两个特征根，第三种情况是当|a[m][m-1]|>E,|a[m-1][m-2]|>E 时，此时需要利用第二种情况得到的矩阵记s=(1)(1)m m a --+()()m m a ，t=(1)(1)m m a --()()m m a -(1)()m m a -(1)()m m a -； 记2 ()()k k k M a sa tI =-+，然后调用QR 分解子程序，求出k K K M Q R =,然后求 1T k K K a Q aQ +=，初始设定第三种情况下进行如上的循环100次，如果无法求出所有的特征 值，则重新设定循环的次数，此程序设定循环次数100次即可满足要求； ③当m=1时，此时要考虑两种情况第一种是若|(1)(0)a | 若|(1)(0)a |>E,此时，由(0)(0)(0)(1) (1)(0)(1)(1) [ ]a a a a 得到的两个特征值即为第二个和第一个特征值； ④当m=0时，此时(0)(0)a 即为第一个特征值； （4）输出QR 分解以后的矩阵； （5）输出矩阵的特征值； 3、求矩阵的特征向量，调用反幂法子程序，求出对应实特征值的特征向量，并输出特征向量。这里反幂算法与上次作业中的相同，子程序亦相似，只是有些变量名稍有改变而已，不再给出具体算法。 二、源程序 #include void nssjh(double temp[10][10]); void sb_qr(double tempa[10][10]); void qr_fangfa(double tempaa[10][10]); void LU(double lujz[10][10]); void FMF(double fmfjz[10][10]); double HK[10][10],MK[10][10],MKM[10][10]; double justin1[10][10],justin2[10][10]; double stupid[10][10]; double fmfu1[10],fmfmiddle[10],fmfy1[10]; double Re[10],Im[10]; //十个特征值的实部和虚部 double malern1,c1; double pig,zdt2,zdt1,sum1,m1,fmmin,fmmin1; double ufs1,ufs11; double u1[10],cp1[10]; double lll[10][10]; double lose[10][10]; double malern; double lul[10][10],luu[10][10]; double cute[10][10]; double fs1,fsqurt1; double llm[10][10]; double forget[10][10]; /*拟上三角化子程序*/ void nssjh(double temp[10][10]) { double c; double u[10],cp[10],e[10]={0}; double fs,fsqurt; for(i=0;i<8;i++) //对于10*10矩阵只需循环8次 { /*将temp矩阵的第i列从i+1个值开始赋给cp[]*/ for(n=0;n<10;n++) { if(n>i) cp[n]=temp[n][i]; else cp[n]=0; } fs=0; //每次循环完了之后重新赋值为0 for(j=0;j<10;j++) { fs+=pow(cp[j],2); } fsqurt=sqrt(fs); if(temp[i+1][i]>0) c=-fsqurt; else c=fsqurt; /*求出u的元素值*/ for(w=0;w<10;w++) { if(w==i+1) u[w]=cp[w]-c; else u[w]=cp[w]; } double ufs=0,ufs1; for(w=0;w<10;w++) { ufs+=pow(u[w],2); } ufs1=0.5*ufs; double ll[10][10]; for(j=0;j<10;j++) { for(n=0;n<10;n++) { if(j==n) ll[j][n]=1-u[j]*u[n]/ufs1; else ll[j][n]=-u[j]*u[n]/ufs1; } } for(j=0;j<10;j++) { for(n=0;n<10;n++) { malern=0; for(w=0;w<10;w++) { malern+=ll[j][w]*temp[w][n]; } if(n==i&&j>i+1&&fabs(malern) lll[j][n]=0; else lll[j][n]=malern; } } double zjbl; for(j=0;j<10;j++) { for(n=0;n<10;n++) { zjbl=0; for(w=0;w<10;w++) { zjbl+=lll[j][w]*ll[w][n]; } temp[j][n]=zjbl; } } } } /*QR分解子程序*/ void qr_fangfa(double tempaa[10][10]) { /*将HK矩阵初值赋为单位矩阵*/ for(i=0;i<10;i++) for(j=0;j<10;j++) { if(i==j) HK[i][j]=1; else HK[i][j]=0; } for(i=0;i { /*将tempaa矩阵的第i列大于等于i行的值赋给cp1[]*/ for(n=0;n { if(n>=i) cp1[n]=tempaa[n][i]; else cp1[n]=0; } fs1=0; //每次循环完了之后重新赋值为0 for(j=0;j { fs1+=pow(cp1[j],2); } fsqurt1=sqrt(fs1); if(tempaa[i][i]>0) c1=-fsqurt1; else c1=fsqurt1; for(w=0;w { if(w==i) u1[w]=cp1[w]-c1; else u1[w]=cp1[w]; } ufs1=0; for(w=0;w { ufs1+=pow(u1[w],2); } ufs11=0.5*ufs1; for(j=0;j { for(n=0;n { if(j==n) llm[j][n]=1-u1[j]*u1[n]/ufs11; else llm[j][n]=-u1[j]*u1[n]/ufs11; } } for(j=0;j for(n=0;n { cute[j][n]=0; //用来储存每次右乘llm矩阵之后的数值 //每次循环开始时初始化为0 } for(j=0;j for(n=0;n { for(w=0;w cute[j][n]+=HK[j][w]*llm[w][n]; } for(j=0;j for(n=0;n { HK[j][n]=cute[j][n]; } for(j=0;j { for(w=0;w { malern1=0; for(n=0;n { malern1+=llm[j][n]*tempaa[n][w]; } if(w==i&&j>i&&fabs(malern1) stupid[j][w]=0; else stupid[j][w]=malern1; } } for(j=0;j for(n=0;n { tempaa[j][n]=stupid[j][n]; } } } /*带双步位移的QR分解*/ void sb_qr(double tempa[10][10]) { double zjl1,zjl2,zjl3; for(m=9;m>=0;) { if(m>1) //第一种情况 { if(fabs(tempa[m][m-1]) { Re[m]=tempa[m][m]; Im[m]=0; m=m-1; } if(fabs(tempa[m][m-1])>E&&fabs(tempa[m-1][m-2]) { zjl1=tempa[m][m]+tempa[m-1][m-1]; zjl2=tempa[m][m]*tempa[m-1][m-1]-tempa[m][m-1]*tempa[m-1][m]; zjl3=pow(zjl1,2)-4*zjl2; if(zjl3>0) { Re[m]=(zjl1+sqrt(zjl3))/2.0; Im[m]=0; Re[m-1]=(zjl1-sqrt(zjl3))/2.0; Im[m-1]=0; } else { Re[m]=0.5*zjl1; Im[m]=sqrt(-zjl3)*0.5; Re[m-1]=0.5*zjl1; Im[m-1]=-sqrt(-zjl3)*0.5; } m=m-2; } if(fabs(tempa[m][m-1])>E&&fabs(tempa[m-1][m-2])>E) //第三种小情况 { int sad; for(sad=0;sad<1000;sad++) { for(i=0;i<10;i++) for(j=0;j<10;j++) { MKM[i][j]=0; //初始化为0矩阵 } zjl1=tempa[m][m]+tempa[m-1][m-1]; zjl2=tempa[m][m]*tempa[m-1][m-1]-tempa[m][m-1]*tempa[m-1][m]; for(i=0;i for(j=0;j { for(w=0;w MKM[i][j]+=tempa[i][w]*tempa[w][j]; } for(i=0;i for(j=0;j { if(i==j) MK[i][j]=MKM[i][j]-zjl1*tempa[i][j]+zjl2; else MK[i][j]=MKM[i][j]-zjl1*tempa[i][j]; } qr_fangfa(MK); //对MK进行qr_fangfa分解,分解完之后的Q矩阵是HK矩阵，R矩阵是MK矩阵 for(i=0;i<10;i++) for(j=0;j<10;j++) { justin1[i][j]=0; //每次都初始化为0 } for(i=0;i<10;i++) for(j=0;j<10;j++) { justin2[i][j]=0; //每次都初始化为0 } for(i=0;i<10;i++) for(j=0;j<10;j++) { justin1[j][i]=HK[i][j]; //作为HK矩阵的转置矩阵 } for(i=0;i<10;i++) for(j=0;j<10;j++) for(w=0;w<10;w++) justin2[i][j]+=justin1[i][w]*tempa[w][j]; for(i=0;i<10;i++) for(j=0;j<10;j++) { lose[i][j]=0; } for(i=0;i<10;i++) for(j=0;j<10;j++) for(w=0;w<10;w++) lose[i][j]+=justin2[i][w]*HK[w][j]; for(i=0;i<10;i++) for(j=0;j<10;j++) { if(i>j+1&&lose[i][j] tempa[i][j]=0; else tempa[i][j]=lose[i][j]; } } } } else if(m==1) //第二种大情况 { if(tempa[1][0] { Re[1]=tempa[1][1]; Im[1]=0; m=0; } else { zjl1=tempa[m][m]+tempa[m-1][m-1]; zjl2=tempa[m][m]*tempa[m-1][m-1]-tempa[m][m-1]*tempa[m-1][m]; zjl3=pow(zjl1,2)-4*zjl2; if(zjl3>0) { Re[1]=(zjl1+sqrt(zjl3))/2.0; Im[1]=0; Re[0]=(zjl1-sqrt(zjl3))/2.0; Im[0]=0; } else { Re[1]=0.5*zjl1; Im[1]=sqrt(-zjl3)*0.5; Re[0]=0.5*zjl1; Im[0]=-sqrt(-zjl3)*0.5; } m=-1; } } else //第三种大情况 { Re[0]=tempa[0][0]; Im[0]=0; m=-1; } } } void LU(double lujz[10][10]) { for(i=0;i<10;i++) for(j=0;j<10;j++) { if(i==j) lul[i][j]=1; else lul[i][j]=0; } for(i=0;i<10;i++) for(j=0;j<10;j++) { luu[i][j]=0; } for(i=0;i<10;i++) { luu[0][i]=lujz[0][i]; //先求出u矩阵的第1行 } for(i=1;i<10;i++) //第1列 { lul[i][0]=lujz[i][0]/luu[0][0]; } for(i=1;i<10;i++) //u第2行 { luu[1][i]=lujz[1][i]-lul[1][0]*luu[0][i]; } for(i=2;i<10;i++) //第2列 { lul[i][1]=(lujz[i][1]-lul[i][0]*luu[0][1])/luu[1][1]; } for(i=2;i<10;i++) //u第3行 { luu[2][i]=lujz[2][i]-lul[2][0]*luu[0][i]-lul[2][1]*luu[1][i]; } for(i=3;i<10;i++) //第3列 { lul[i][2]=(lujz[i][2]-lul[i][0]*luu[0][2]-lul[i][1]*luu[1][2])/luu[2][2]; } for(i=3;i<10;i++) //u第4行 { luu[3][i]=lujz[3][i]-lul[3][0]*luu[0][i]-lul[3][1]*luu[1][i]-lul[3][2]*luu[2][i]; } for(i=4;i<10;i++) //第4列 { lul[i][3]=(lujz[i][3]-lul[i][0]*luu[0][3]-lul[i][1]*luu[1][3] -lul[i][2]*luu[2][3])/luu[3][3]; } for(i=4;i<10;i++) //u第5行 { luu[4][i]=lujz[4][i]-lul[4][0]*luu[0][i]-lul[4][1]*luu[1][i] -lul[4][2]*luu[2][i]-lul[4][3]*luu[3][i]; } for(i=5;i<10;i++) //第5列 { lul[i][4]=(lujz[i][4]-lul[i][0]*luu[0][4]-lul[i][1]*luu[1][4] -lul[i][2]*luu[2][4]-lul[i][3]*luu[3][4])/luu[4][4]; } for(i=5;i<10;i++) //u第6行 { luu[5][i]=lujz[5][i]-lul[5][0]*luu[0][i]-lul[5][1]*luu[1][i] -lul[5][2]*luu[2][i]-lul[5][3]*luu[3][i]-lul[5][4]*luu[4][i]; } for(i=6;i<10;i++) //第6列 { lul[i][5]=(lujz[i][5]-lul[i][0]*luu[0][5]-lul[i][1]*luu[1][5] -lul[i][2]*luu[2][5]-lul[i][3]*luu[3][5]- lul[i][4]*luu[4][5])/luu[5][5]; } for(i=6;i<10;i++) //u第7行 { luu[6][i]=lujz[6][i]-lul[6][0]*luu[0][i]-lul[6][1]*luu[1][i] -lul[6][2]*luu[2][i]-lul[6][3]*luu[3][i]- lul[6][4]*luu[4][i]-lul[6][5]*luu[5][i]; } for(i=7;i<10;i++) //第7列 { lul[i][6]=(lujz[i][6]-lul[i][0]*luu[0][6]-lul[i][1]*luu[1][6] -lul[i][2]*luu[2][6]-lul[i][3]*luu[3][6]- lul[i][4]*luu[4][6]-lul[i][5]*luu[5][6])/luu[6][6]; } for(i=7;i<10;i++) //u第8行 luu[7][i]=lujz[7][i]-lul[7][0]*luu[0][i]-lul[7][1]*luu[1][i] -lul[7][2]*luu[2][i]-lul[7][3]*luu[3][i]- lul[7][4]*luu[4][i]-lul[7][5]*luu[5][i] -lul[7][6]*luu[6][i]; } for(i=8;i<10;i++) //第8列 { lul[i][7]=(lujz[i][7]-lul[i][0]*luu[0][7]-lul[i][1]*luu[1][7] -lul[i][2]*luu[2][7]-lul[i][3]*luu[3][7]- lul[i][4]*luu[4][7]-lul[i][5]*luu[5][7]- lul[i][6]*luu[6][7])/luu[7][7]; } for(i=8;i<10;i++) //u第9行 { luu[8][i]=lujz[8][i]-lul[8][0]*luu[0][i]-lul[8][1]*luu[1][i] -lul[8][2]*luu[2][i]-lul[8][3]*luu[3][i]- lul[8][4]*luu[4][i]-lul[8][5]*luu[5][i] -lul[8][6]*luu[6][i]-lul[8][7]*luu[7][i]; } for(i=9;i<10;i++) //第9列 { lul[i][8]=(lujz[i][8]-lul[i][0]*luu[0][8]-lul[i][1]*luu[1][8] -lul[i][2]*luu[2][8]-lul[i][3]*luu[3][8]- lul[i][4]*luu[4][8]-lul[i][5]*luu[5][8]- lul[i][6]*luu[6][8]-lul[i][7]*luu[7][8])/luu[8][8]; } luu[9][9]=lujz[9][9]-lul[9][0]*luu[0][9]-lul[9][1]*luu[1][9] //u第10行 -lul[9][2]*luu[2][9]-lul[9][3]*luu[3][9]- lul[9][4]*luu[4][9]-lul[9][5]*luu[5][9] -lul[9][6]*luu[6][9]-lul[9][7]*luu[7][9]- lul[9][8]*luu[8][9]; } /*反幂法子程序：用来求解矩阵按模最小的特征值和特征向量*/ void FMF(double fmfjz[10][10]) { /*u1矩阵赋初值,初始值都是1*/ for(i=0;i<10;i++) { fmfu1[i]=1; zdt2=0; //先假设zdt2=0 do //在未达到精度时进行循环 { zdt1=zdt2; //将一轮计算后的结果赋给zdt1 sum1=0; for(i=0;i<10;i++) //u矩阵元素求平方和 { sum1+=pow(fmfu1[i],2); } m1=sqrt(sum1); for(i=0;i<10;i++) { fmfy1[i]=fmfu1[i]/m1; //y(k-1) } fmfmiddle[0]=fmfy1[0]; fmfmiddle[1]=fmfy1[1]-lul[1][0]*fmfmiddle[0]; fmfmiddle[2]=fmfy1[2]-lul[2][0]*fmfmiddle[0]-lul[2][1]*fmfmiddle[1]; fmfmiddle[3]=fmfy1[3]-lul[3][0]*fmfmiddle[0]-lul[3][1]*fmfmiddle[1] -lul[3][2]*fmfmiddle[2]; fmfmiddle[4]=fmfy1[4]-lul[4][0]*fmfmiddle[0]-lul[4][1]*fmfmiddle[1] -lul[4][2]*fmfmiddle[2]-lul[4][3]*fmfmiddle[3]; fmfmiddle[5]=fmfy1[5]-lul[5][0]*fmfmiddle[0]-lul[5][1]*fmfmiddle[1] -lul[5][2]*fmfmiddle[2]-lul[5][3]*fmfmiddle[3]- lul[5][4]*fmfmiddle[4]; fmfmiddle[6]=fmfy1[6]-lul[6][0]*fmfmiddle[0]-lul[6][1]*fmfmiddle[1] -lul[6][2]*fmfmiddle[2]-lul[6][3]*fmfmiddle[3]- lul[6][4]*fmfmiddle[4]-lul[6][5]*fmfmiddle[5]; fmfmiddle[7]=fmfy1[7]-lul[7][0]*fmfmiddle[0]-lul[7][1]*fmfmiddle[1] -lul[7][2]*fmfmiddle[2]-lul[7][3]*fmfmiddle[3]- lul[7][4]*fmfmiddle[4]-lul[7][5]*fmfmiddle[5]- lul[7][6]*fmfmiddle[6]; fmfmiddle[8]=fmfy1[8]-lul[8][0]*fmfmiddle[0]-lul[8][1]*fmfmiddle[1] -lul[8][2]*fmfmiddle[2]-lul[8][3]*fmfmiddle[3]- lul[8][4]*fmfmiddle[4]-lul[8][5]*fmfmiddle[5]- lul[8][6]*fmfmiddle[6]-lul[8][7]*fmfmiddle[7]; fmfmiddle[9]=fmfy1[9]-lul[9][0]*fmfmiddle[0]-lul[9][1]*fmfmiddle[1] -lul[9][2]*fmfmiddle[2]-lul[9][3]*fmfmiddle[3]- lul[9][4]*fmfmiddle[4]-lul[9][5]*fmfmiddle[5]- lul[9][6]*fmfmiddle[6]-lul[9][7]*fmfmiddle[7]- lul[9][8]*fmfmiddle[8]; fmfu1[9]=fmfmiddle[9]; fmfu1[8]=(fmfmiddle[8]-luu[8][9]*fmfu1[9])/luu[8][8]; fmfu1[7]=(fmfmiddle[7]-luu[7][9]*fmfu1[9]-luu[7][8]*fmfu1[8])/luu[7][7]; fmfu1[6]=(fmfmiddle[6]-luu[6][9]*fmfu1[9]-luu[6][8]*fmfu1[8] -luu[6][7]*fmfu1[7])/luu[6][6]; fmfu1[5]=(fmfmiddle[5]-luu[5][9]*fmfu1[9]-luu[5][8]*fmfu1[8] -luu[5][7]*fmfu1[7]-luu[5][6]*fmfu1[6])/luu[5][5]; fmfu1[4]=(fmfmiddle[4]-luu[4][9]*fmfu1[9]-luu[4][8]*fmfu1[8] -luu[4][7]*fmfu1[7]-luu[4][6]*fmfu1[6] -luu[4][5]*fmfu1[5])/luu[4][4]; fmfu1[3]=(fmfmiddle[3]-luu[3][9]*fmfu1[9]-luu[3][8]*fmfu1[8] -luu[3][7]*fmfu1[7]-luu[3][6]*fmfu1[6] -luu[3][5]*fmfu1[5]-luu[3][4]*fmfu1[4])/luu[3][3]; fmfu1[2]=(fmfmiddle[2]-luu[2][9]*fmfu1[9]-luu[2][8]*fmfu1[8] -luu[2][7]*fmfu1[7]-luu[2][6]*fmfu1[6] -luu[2][5]*fmfu1[5]-luu[2][4]*fmfu1[4]- luu[2][3]*fmfu1[3])/luu[2][2]; fmfu1[1]=(fmfmiddle[1]-luu[1][9]*fmfu1[9]-luu[1][8]*fmfu1[8] -luu[1][7]*fmfu1[7]-luu[1][6]*fmfu1[6] -luu[1][5]*fmfu1[5]-luu[1][4]*fmfu1[4]- luu[1][3]*fmfu1[3]-luu[1][2]*fmfu1[2])/luu[1][1]; fmfu1[0]=(fmfmiddle[0]-luu[0][9]*fmfu1[9]-luu[0][8]*fmfu1[8] -luu[0][7]*fmfu1[7]-luu[0][6]*fmfu1[6] -luu[0][5]*fmfu1[5]-luu[0][4]*fmfu1[4]- luu[0][3]*fmfu1[3]-luu[0][2]*fmfu1[2]- luu[0][1]*fmfu1[1])/luu[0][0]; double fmsum=0; for(i=0;i<10;i++) { fmsum+=fmfu1[i]*fmfy1[i]; } zdt2=fmsum; }while(fabs(zdt2-zdt1)/fabs(zdt2)>E); } /*主程序主要用来求特征值和对应于实特征值的特征向量*/ void main() { /*初始化矩阵*/ for(i=0;i<10;i++) for(j=0;j<10;j++) { if(i==j) a[i][j]=1.5*cos((i+1)+1.2*(j+1)); else a[i][j]=sin(0.5*(i+1)+0.2*(j+1)); } cout< cout.precision(11); cout.unsetf(ios::fixed); nssjh(a); //调用拟上三角化矩阵 /*输出拟上三角矩阵*/ for(j=0;j<10;j++) for(n=0;n<10;n++) { cout< if(n!=0&&n%9==0) cout< } cout< sb_qr(a); //调用带双步位移的qr_fangfa分解子程序求特征值 for(j=0;j<10;j++) for(n=0;n<10;n++) { cout< if(n!=0&&n%9==0) cout< } cout< cout<<"Re[0]="< < < < < < < < < < < for(i=0;i<10;i++) //将a矩阵复原 for(j=0;j<10;j++) { if(i==j) a[i][j]=1.5*cos((i+1)+1.2*(j+1)); else a[i][j]=sin(0.5*(i+1)+0.2*(j+1)); } for(j=0;j<10;j++) for(w=0;w<10;w++) { if(j==w) a[j][w]=a[j][w]-Re[3]; else a[j][w]=a[j][w]; } LU(a); FMF(a); cout<<"与特征值"< < < < < < < < < < < for(i=0;i<10;i++) //将a矩阵复原 for(j=0;j<10;j++) { if(i==j) a[i][j]=1.5*cos((i+1)+1.2*(j+1)); else a[i][j]=sin(0.5*(i+1)+0.2*(j+1)); } for(j=0;j<10;j++) for(w=0;w<10;w++) { if(j==w) a[j][w]=a[j][w]-Re[6]; else 数值分析 第二章 2．当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) ()()3 x 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 ==-===-=--==-+-----==------= =-+-- 则二次拉格朗日插值多项式为 2 20 ()()k k k L x y l x ==∑ 0223()4() 14 (1)(2)(1)(1)23 537623 l x l x x x x x x x =-+=---+ -+= +- 6．设,0,1,,j x j n =L 为互异节点，求证： （1） 0()n k k j j j x l x x =≡∑ (0,1,,);k n =L （2）0 ()()0n k j j j x x l x =-≡∑ (0,1,,);k n =L 证明 （1） 令()k f x x = 若插值节点为,0,1,,j x j n =L ，则函数()f x 的n 次插值多项式为0 ()()n k n j j j L x x l x == ∑。 插值余项为(1)1() ()()()()(1)! n n n n f R x f x L x x n ξω++=-= + 又,k n ≤Q (1)()0 ()0 n n f R x ξ+∴=∴= 0()n k k j j j x l x x =∴=∑ (0,1,,);k n =L 0 000 (2)()() (())()()(()) n k j j j n n j i k i k j j j i n n i k i i k j j i j x x l x C x x l x C x x l x =-==-==-=-=-∑∑∑∑∑ 0i n ≤≤Q 又 由上题结论可知 ()n k i j j j x l x x ==∑ ()()0 n i k i i k i k C x x x x -=∴=-=-=∑原式 ∴得证。 7设[]2 (),f x C a b ∈且()()0,f a f b ==求证： 21 max ()()max ().8 a x b a x b f x b a f x ≤≤≤≤''≤- 解：令01,x a x b ==，以此为插值节点，则线性插值多项式为 10 101010 ()() ()x x x x L x f x f x x x x x --=+-- =() () x b x a f a f b a b x a --=+-- 1()()0()0 f a f b L x ==∴=Q 又 插值余项为1011 ()()()()()()2 R x f x L x f x x x x x ''=-= -- 011 ()()()()2 f x f x x x x x ''∴= -- 《数值分析B》大作业一 SY1103120 朱舜杰 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求 出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和 最小特征值。 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。 二．源程序 #include 第一章 1、ln2＝0.69314718…，精确到 10－3 的近似值是多少？ 解 精确到 10－3＝0.001，即绝对误差限是 e ＝0.05％，故至少要保留小数点后三位才可以。 ln2≈0.693。 2、设115.80,1025.621≈≈x x 均具有5位有效数字，试估计由这些数据计算21x x ， 21x x +的绝对误差限 解：记126.1025, 80.115x x == 则有11232411 10, | 102|||2 x x x x --≤?-≤?- 所以 121212121212211122||||||||||||x x x x x x x x x x x x x x x x x x -=-+-+≤-- 3411 80.11610 6.10102522 0.007057-==??+≤?? 1212112243|()|||11 |10100.0005522 |x x x x x x x x --≤≤?+?=+-+-+- 3、一个园柱体的工件,直径d 为10.250.25mm,高h 为40.00 1.00mm,则它的体 积V 的近似值、误差和相对误差为多少。 解： ()() 22222222 4 314210254000000330064 221025400002510251002436444 3300624362436 0073873833006 , .....; ()()()......, ..().()..% .r d h V d h V mm d h V dh d d h V mm V V V πππππεεεεε= ≈=??===+=???+?==±====第二章： 1、分别利用下面四个点的Lagrange 插值多项式和Newton 插值多项式N 3(x )， 计算L 3(0.5)及N 3(-0.5) x －2 －1 0 1 f (x ) －1 1 2 数值分析作业答案 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange插值基底。 （3）用Newton基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为: ， 所以： 所以f(x)的二次插值多项式为： （2）用Lagrange插值基底 Lagrange插值多项式为： 所以f(x)的二次插值多项式为： (3) 用Newton基底: 均差表如下： xk f(xk) 一阶均差二阶均差 1 0 -1 -3 3/2 2 4 7/ 3 5/6 Newton插值多项式为： 所以f(x)的二次插值多项式为： 由以上计算可知，三种方法得到的多项式是相同的。 6、在上给出的等距节点函数表，若用二次插值求ex的近似值，要使截断误差不超过10-6，问使用函数表的步长h应取多少？ 解：以xi-1,xi,xi+1为插值节点多项式的截断误差，则有 式中 令得 插值点个数 是奇数，故实际可采用的函数值表步长 8、，求及。 解：由均差的性质可知，均差与导数有如下关系： 所以有： 15、证明两点三次Hermite插值余项是 并由此求出分段三次Hermite插值的误差限。 证明：利用[xk，xk+1]上两点三次Hermite插值条件 知有二重零点xk和k+1。设 确定函数k(x): 当或xk+1时k(x)取任何有限值均可； 当时，，构造关于变量t的函数 显然有 在[xk,x][x,xk+1]上对g(x)使用Rolle定理，存在及使得 在，，上对使用Rolle定理，存在，和使得 再依次对和使用Rolle定理，知至少存在使得 而，将代入，得到 推导过程表明依赖于及x 综合以上过程有： 确定误差限： 记为f(x)在[a,b]上基于等距节点的分段三次Hermite插值函数。在区间[xk，xk+1]上有 而最值 进而得误差估计： 16、求一个次数不高于4次的多项式，使它满足，，。 第二章复习与思考题 1.什么是拉格朗日插值基函数？它们是如何构造的？有何重要性质？ 答：若n 次多项式()),,1,0(n j x l j =在1+n 个节点n x x x <<< 10上满足条件 (),,,1,0,, ,0, ,1n k j j k j k x l k j =?? ?≠== 则称这1+n 个n 次多项式()()()x l x l x l n ,,,10 为节点n x x x ,,,10 上的n 次拉格朗日插值基函数. 以()x l k 为例，由()x l k 所满足的条件以及()x l k 为n 次多项式，可设 ()()()()()n k k k x x x x x x x x A x l ----=+- 110， 其中A 为常数，利用()1=k k x l 得 ()()()()n k k k k k k x x x x x x x x A ----=+- 1101， 故 ()()()() n k k k k k k x x x x x x x x A ----= +- 1101 ， 即 ()()()()()()()()∏ ≠=+-+---=--------=n k j j j k j n k k k k k k n k k k x x x x x x x x x x x x x x x x x x x x x l 0110110)( . 对于()),,1,0(n i x l i =，有 ()n k x x l x n i k i k i ,,1,00 ==∑=，特别当0=k 时，有 ()∑==n i i x l 0 1. 2.什么是牛顿基函数？它与单项式基{ }n x x ,,,1 有何不同？ 答：称()()()(){ }10100,,,,1------n x x x x x x x x x x 为节点n x x x ,,,10 上的牛顿基函数，利用牛顿基函数，节点n x x x ,,,10 上的n 次牛顿插值多项式()x P n 可以表示为 ()()()()10010---++-+=n n n x x x x a x x a a x P 其中[]n k x x x f a k k ,,1,0,,,,10 ==.与拉格朗日插值多项式不同，牛顿插值基函数在增加节点时可以通过递推逐步得到高次的插值多项式，例如 ()()()()k k k k x x x x a x P x P --+=++ 011, 数值分析习题集 （适合课程《数值方法A 》和《数值方法B 》） 长沙理工大学 第一章 绪 论 1. 设x >0,x 的相对误差为δ,求ln x 的误差. 2. 设x 的相对误差为2％,求n x 的相对误差. 3. 下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指出 它们是几位有效数字: *****123451.1021,0.031,385.6,56.430,7 1.0.x x x x x =====? 4. 利用公式(3.3)求下列各近似值的误差限: ********12412324(),(),()/,i x x x ii x x x iii x x ++其中**** 1234 ,,,x x x x 均为第3题所给的数. 5. 计算球体积要使相对误差限为1％,问度量半径R 时允许的相对误差限是多少? 6. 设028,Y =按递推公式 1n n Y Y -=( n=1,2,…) 计算到100Y .27.982(五位有效数字),试问计算100Y 将有多大误差? 7. 求方程2 5610x x -+=的两个根,使它至少具有四位有效数字27.982). 8. 当N 充分大时,怎样求2 1 1N dx x +∞+?? 9. 正方形的边长大约为100㎝,应怎样测量才能使其面积误差不超过1㎝2 ? 10. 设 212S gt = 假定g 是准确的,而对t 的测量有±0.1秒的误差,证明当t 增加时S 的绝对 误差增加,而相对误差却减小. 11. 序列 {}n y 满足递推关系1101n n y y -=-(n=1,2,…),若0 1.41y =≈(三位有效数字), 计算到 10y 时误差有多大?这个计算过程稳定吗? 12. 计算6 1)f =, 1.4≈,利用下列等式计算,哪一个得到的结果最好? 3 -- 13. ()ln(f x x =,求f (30)的值.若开平方用六位函数表,问求对数时误差有多大?若 目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。 这样即求出所有A所有实特征值对应的一个特征向量。 #include ￥ 数值分析思考题1 1、讨论绝对误差（限）、相对误差（限）与有效数字之间的关系。 2、相对误差在什么情况下可以用下式代替 3、查阅何谓问题的“病态性”，并区分与“数值稳定性”的不同点。 4、取 ，计算 ，下列方法中哪种最好为什么（1）(3 3-，（2）(2 7-，（3） ()3 1 3+ ，（4） ()6 1 1 ，（5）99- , 数值实验 数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。 Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。 数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 $ r e x x e x x ** * ** - == 141 . ≈)61 《数值分析》 重点考察内容及各章作业答案 学院： 学号： 姓名： 重点考察内容 基本概念（收敛阶，收敛条件，收敛区域等）, 简单欧拉法。 第一章基础 掌握：误差的种类，截断误差，舍入误差的来源，有效数字的判断。 了解：误差限，算法及要注意的问题。 第二章插值 掌握：Hermite插值，牛顿插值，差商计算，插值误差估计。 了解：Lagrange插值 第三章数据拟合 掌握：给出几个点求线性拟合曲线。 了解：最小二乘原理 第四章数值积分微分 掌握：梯形公式，Simpson公式，代数精度，Gauss积分，带权Gauss积分公式推导，复化梯形公式推导及算法。 了解：数值微分，积分余项 第五章直接法 掌握：LU分解求线性方程组，运算量 了解：Gauss消去法，LDL，追赶法 第六章迭代法 掌握：Jacobi,Gauss-Seidel迭代格式构造，敛散性分析，向量、矩阵的范数、谱半径 了解：SOR迭代 第七章Nolinear迭代法 掌握：牛顿迭代格式构造，简单迭代法构造、敛散性分析，收敛阶。 了解：二分法，弦截法 第八章ODE解法 掌握：Euler公式构造、收敛阶。 了解：梯形Euler公式、收敛阶，改进Euler公式 题目类型：填空，计算，证明综合题 第一章 误差 1. 科学计算中的误差来源有4个，分别是________，________，________，________。 2. 用Taylor 展开近似计算函数000()()'()()f x f x f x x x ≈+-，这里产生是什么误差？ 3. 0.7499作 3 4 的近似值，是______位有效数字，65.380是舍入得到的近似值，有____几位有效数字，相对误差限为_______. 0.0032581是四舍五入得到的近似值，有_______位有效数字. 4. 改变下列表达式，使计算结果比较精确： （1）11,||1121x x x x --++ （2 ||1x （3） 1cos ,0,|| 1.x x x x -≠ （4）sin sin ,αβαβ-≈ 5. 采用下列各式计算61)时，哪个计算效果最好？并说明理由。 （1） （2 ）99-3 ）6 (3-（4 6. 已知近似数*x 有4位有效数字，求其相对误差限。 上机实验题： 1、利用Taylor 展开公式计算0! k x k x e k ∞ ==∑，编一段小程序，上机用单精度计算x e 的函数 值. 分别取x =1，5，10，20，-1，-5，-10，-15，-20，观察所得结果是否合理，如不合理请分析原因并给出解决方法. 2、已知定积分1 ,0,1,2,,206 n n x I dx n x ==+? ,有如下的递推关系 111 110 0(6)61666 n n n n n x x x x I dx dx I x x n ---+-===++-? ? 可建立两种等价的计算公式 (1) 1016,0.154n n I I I n -= -=取； (2) 12011),0.6n n I nI I n -=-=（取 数值分析大作业（2013年5月） 金洋洋（12721512），机自系 1.下列各数都是经过四舍五入得到的近似值，试分别指出它 们的绝对误差限, 相对误差限和有效数字的位数。 X1 =5.420, x 2 =0.5420, x 3=0.00542, x 4 =6000, x 5=50.610? 解：根据定义：如果*x 的绝对误差限 不超过x 的某个数位的半个单位，则从*x 的首位非零数字到该位都是有效数字。 显然根据四舍五入原则得到的近视值，全部都是有效数字。 因而在这里有：n1=4, n2=4, n3=3, n4=4, n5=1 （n 表示x 有效数字的位数） 对x1：有a1=5, m1=1 （其中a1表示x 的首位非零数字，m1表示x1的整数位数） 所以有绝对误差限 143 11 (1)101022 x ε--≤ ?=? 相对误差限 31() 0.510(1)0.00923%5.4201 r x x x εε-?= == 对x2：有a2=5, m2=0 所以有绝对误差限 044 11 (2)101022 x ε--≤ ?=? 相对误差限 42() 0.510(2)0.00923%0.54202 r x x x εε-?= == 对x3：有a3=5, m3=-2 所以有绝对误差限 235 11 (3)101022 x ε---≤ ?=? 相对误差限 53() 0.510(3)0.0923%0.005423 r x x x εε-?= == 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022 x ε-≤?= 相对误差限 4() 0.5 (4)0.0083%6000 4 r x x x εε= = = 对x5：有a5=6, m5=5 所以有绝对误差限 514 11(5)101022 x ε-≤ ?=? 相对误差限 45() 0.510(5)8.3%600005 r x x x εε?= == 1、讨论绝对误差（限）、相对误差（限）与有效数字之间的关系。 2、相对误差在什么情况下可以用下式代替？ 3、查阅何谓问题的“病态性”，并区分与“数值稳定性”的不同点。 4、 取 ，计算 ，不用计算而直接判断下列式子中哪 种计算效果最好？为什么？ （1）(3 3-，（2）(2 7-，（3） (3 1 3+，（4） ) 6 11 ，（5）99-5. 应用梯形公式 ))()((2b f a f a b T +-= 计算积分1 0x I e dx -=?的近似值，在整个计算过程中按四舍五入规则取五位小数。计算中产生的误差的主要原因是截断误差还是舍入误差？为什么？ 6. 下列各数都是经过四舍五入得到的近似值，试指出他们有几位有效数字，并给出其绝对误差限与相对误差限。 (1) 1021.1*1=x ；(2) 031.0*2=x ；(3) 40.560*3=x 。 7. 下列公式如何计算才比较准确？ (1) 212 x e -，1x <<；(2) 12 1 N N dx x ++? ，1>>N ；(3) ，1x >>。 8. 序列{}n y 满足递推关系1101n n y y -=-，12,,n =，若0141.y =≈，计算到10y 时误差有多大？这个计算过程数值稳定吗？ r e x x e x x ***** -== 141.≈) 6 1 1、怎样确定一个隔根区间？如何求解一个方程的全部实根？如：已知方程：1020()x f x e x =+-=在(),-∞+∞有实数根，用二分法求它的全部实根，要求误差满足210*k x x --<？若要求6*10k x x --<，需二分区 间多少次？ 2、求解一个非线性方程的迭代法有哪些充分条件可以保障迭代序列收敛于方程的根？对方程3210()f x x x =--=，试构造两种不同的迭代法，且均收敛于方程在[]12,中的唯一根。 3、设0a >，应用牛顿法于方程30x a -= 确定常数,p q 和r 使得迭代法 2 125k k k k qa ra x px x x +=++， 012,, , k = 4、对于不动点方程()x x ?=，()x ?满足映内性和压缩性是存在不动点的充分条件，他们也是必要条件吗？试证明：（1）函数21()x x ?=-在闭区间[]02,上不是映内的，但在其上有不动点；（2）函数 1()ln()x x e ?=+在任何区间[],a b 上都是压缩的，但没有不动点。 5、设*x 是方程0()f x =的根，且0*'()f x ≠，''()f x 在*x 的某个邻域上连续。试证明：Newton 迭代序列{}k x 满足 12122**()''() lim () '()k k k k k x x f x x x f x -→∞---=-- 6. 设有方程1 12 sin x x =+。对于迭代法1112 ()sin()k k k x x x ?+==+，试证：对 任何15.b ≥，迭代函数()x ?在闭区间[0.5,b]上满足映内性和压缩性。用所给方 习题1 1． 填空题 (1) 为便于算法在计算机上实现，必须将一个数学问题分解为 的 运算； (2) 在数值计算中为避免损失有效数字，尽量避免两个 数作减法运算；为避免 误差的扩大，也尽量避免分母的绝对值 分子的绝对值； (3) 误差有四大来源,数值分析主要处理其中的 和 ； (4) 有效数字越多,相对误差越 ； 2. 用例1.4的算法计算10,迭代3次,计算结果保留4位有效数字. 3. 推导开平方运算的误差限公式,并说明什么情况下结果误差不大于自变量误差. 4. 以下各数都是对准确值进行四舍五入得到的近似数,指出它们的有效数位、误差限和相对误差限. 95123450304051104000003346087510., ., , ., .x x x x x -==?===? 5. 证明1.2.3之定理1.1. 6. 若钢珠的的直径d 的相对误差为1.0%,则它的体积V 的相对误差将为多少。（假定钢珠为标准的球形） 7. 若跑道长的测量有0.1%的误差,对400m 成绩为60s 的运动员的成绩将会带来多大的误差和相对误差. 8. 为使20的近似数相对误差小于0.05%,试问该保留几位有效数字. 9. 一个园柱体的工件,直径d 为10.25±0.25mm,高h 为40.00±1.00mm,则它的体积V 的近似值、误差和相对误差为多少． 10 证明对一元函数运算有 r r xf x f x k x k f x εε'≈= () (())(),() 其中 并求出157f x x x ==()tan ,.时的k 值，从而说明f x x =()tan 在2 x π ≈时是病态问题． 11. 定义多元函数运算 1 1 1,,(),n n i i i i i i S c x c x εε====≤∑∑其中 求出S ε()的表达式，并说明i c 全为正数时，计算是稳定的，i c 有正有负时，误差难以控制． 12. 下列各式应如何改进,使计算更准确： 《数值分析》计算实习报告 第二题 院系：机械工程及自动化学院 学号： 姓名： 2017年11月 一、题目要求 试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知 a ij ={ sin (0.5i +0.2j ) i ≠j 1.52cos (i +1.2j ) i =j (i,j =1,2, (10) 说明： 1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10?12。 2.打印以下内容： （1）全部源程序； （2）矩阵A 经过拟上三角化后所得的矩阵A (n?1)； （3）对矩阵A (n?1)实行QR 方法迭代结束后所得的矩阵； （4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,?,10)，其中R i =Re(λi )，I i = Im(λi ) 。若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。 3.采用e 型数输出实型数，并且至少显示12位有效数字。 二、算法设计思路和方案 1. 将矩阵A 拟上三角化得到矩阵A (n?1) 为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n?1)，然后用QR 方法计算A (n?1)的全部特征值，而A (n?1)的特征值就是A 的特征值。具体算法如下： 记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2, ,;,1,,)ij a i n j r r n ==+。 对于1,2,,2r n =-执行 （1）若() (2,3,,)r ir a i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。 第二章思考题 1.画出残差校正法的程序框图 2.查阅目前解决病态线性方程组的方法有哪些，并计较其优劣。 答：对病态线性方程组的求解可采用一下方法求解： （1）采用高精度的算术运算。采用双精度，可改善和减轻病态矩阵的影响。 （2）平衡方法。当A 的元素的数量级差别很大时，采用行平衡或列平衡的方法可降低A 的条件数。 （3）残差矫正法。设A 非奇异且Cond （A ）不特别大，方程组Ax=b 病态但不特别严重，这时可用残差矫正法求解Ax=b 。 （4）残量迭代法。残量迭代法是在用Gauss 消去法求出方程组AX b =的近似解后，进行以下的计算： PA LU =. 然后，重复迭代： b AX γ=- LY P γ= UZ Y = X X Z =+ 其中，PA LU =，X ，Y 和Z 用t 位有效数字计算，b AX -用2t 位有效数字生成。 残量迭代法的计算量比较少，在机器上非常容易执行，而且计算精度也比较高，但是它不是对所有的病态方程组都适用，当()10t Cond A ≥时，其计算结果就不够准确。 （5）加权迭代改善法。 加权迭代改善法是对方程组AX b =构造一个迭代过程：()()()1k k A I X b X λλ++=+ ，其中0λ≠为一常数，I 为与A 同阶的单位方阵,0,1,2,k = 为迭代次数，()0X 为解的初始值，()k X 是第k 次迭代后求得的近似解。只要λ取得合适，A I λ+的逆矩阵()1 A I λ-+便存在。 加权迭代改善法不必选主元且保持原系数矩阵的稀疏性，通过加权因子λ的选取来改善矩阵的条件数，得到了比较好的计算结果，但是由于加权因子λ的选取也使得加权迭代改善法的应用受到了限制。 （6）误差转移法。误差转移法是基于即使方程组的计算解的精度不高，也可获得相对较小的余量这一特点而设计的。 设方程组AX b =的计算解为x τ，既然b Ax τ=对误差很大的解x τ也能比较准确的成立，因而，如果求解r x cy *=其中，c 为n n ?阶非奇异矩阵，则即使计算解r y 的误差比较大，得到比较准确的x *。， 在这种解法中，问题的病态性固然会导致解的巨大误差，但这种误差直接反映在r y 上，对x *的影响则小得多，因为主要的误差已经从原来的x *转移到中间量r y 上了。 误差转移法的原理及实现都十分简捷，仅运用了常规的行列均衡 数值分析第二次大作业 史立峰 SY1505327 一、 方案 （1）利用循环结构将sin(0.50.2)() 1.5cos( 1.2)() {i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的 矩阵A ； （2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。 对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij ,,1,;,,2,1) ( +==。 对于2,,2,1-=n r 执行 1. 若 ()n r r i a r ir ,,3,2) ( ++=全为零，则令A(r+1) =A(r),转5；否则转2。 2. 计算 () ∑+== n r i r ir r a d 1 2 )( ()( )r r r r r r r r r r d c a d a c ==-=++则取,0sgn ) (,1)(,1若 )(,12r r r r r r a c c h +-= 3. 令 () n T r nr r r r r r r r r R a a c a u ∈-=++) ()(,2)(,1,,,,0,,0 。 4. 计算 r r T r r h u A p /)(= r r r r h u A q /)(= r r T r r h u p t /= r r r r u t q -=ω T r r T r r r r p u u A A --=+ω)()1( 5. 继续。 （3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下： 1. 给定精度水平0>ε和迭代最大次数L 。 2. 记n n ij n a A A ?-==][) 1()1()1(，令n m k ==,1。 第一章 绪论 姓名 学号 班级 习题主要考察点：有效数字的计算、计算方法的比较选择、误差和误差限的计算。 1 若误差限为5105.0-?,那么近似数有几位有效数字（有效数字的计算） 解：2*103400.0-?=x ，325*102 1102 1---?=?≤-x x 故具有3位有效数字。 2 14159.3=π具有4位有效数字的近似值是多少（有效数字的计算） 解：10314159.0?= π，欲使其近似值*π具有4位有效数字，必需 41*102 1 -?≤-ππ，3*3102 1102 1--?+≤≤?-πππ，即14209.314109.3*≤≤π 3 已知2031.1=a ，978.0=b 是经过四舍五入后得到的近似值，问b a +， b a ?有几位有效数字（有效数字的计算） 解：3*1021 -?≤-a a ，2*102 1-?≤-b b ，而1811.2=+b a ，1766.1=?b a 2123****102 1 10211021)()(---?≤?+?≤ -+-≤+-+b b a a b a b a 故b a +至少具有2位有效数字。 2123*****102 1 0065.01022031.1102978.0)()(---?≤=?+?≤ -+-≤-b b a a a b b a ab 故b a ?至少具有2位有效数字。 4 设0>x ，x 的相对误差为δ，求x ln 的误差和相对误差（误差的计算） 解：已知δ=-* *x x x ，则误差为 δ=-= -* **ln ln x x x x x 则相对误差为 * * ** * * ln ln 1ln ln ln x x x x x x x x δ = -= - 5测得某圆柱体高度h 的值为cm h 20*=，底面半径r 的值为cm r 5*=， 已知cm h h 2.0||*≤-，cm r r 1.0||*≤-，求圆柱体体积h r v 2π=的绝对误差 限与相对误差限。（误差限的计算） 解：*2******2),(),(h h r r r h r r h v r h v -+-≤-ππ 绝对误差限为 πππ252.051.02052)5,20(),(2=??+????≤-v r h v 相对误差限为 %420 1 20525) 5,20() 5,20(),(2 ==??≤ -ππv v r h v 6 设x 的相对误差为%a ,求n x y =的相对误差。（函数误差的计算） 解：%* *a x x x =-， )%(* **** *na x x x n x x x y y y n n n =-≤-= - 7计算球的体积，为了使体积的相对误差限为%1，问度量半径r 时允许的相对误差限为多大（函数误差的计算） 第2章线性方程组的解法 --------学习小结 一、本章学习体会 通过本章知识的学习我首先了解到求解线性方程组的方法可分为两类：直接法和迭代法。计算机虽然运行速度很快，但面对运算量超级多的问题，计算机还是需要很长的时间进行运算，所以，确定快捷精确的求解线性方程组的方法是非常必要的。 本章分为四个小节，其中前两节Gauss消去法和直接三角分解法因为由之前《线性代数》学习的一定功底，学习起来还较为简单，加之王老师可是的讲解与习题测试，对这一部分有了较好的掌握。第三节矩阵的条件数与病态方程组，我 Ax 的系数矩阵A与左端向量b的元素往往是通首先了解到的是线性方程组b 过观测或计算而得到，因而会带有误差。即使原始数据是精确的，但存放到计算机后由于受字长的限制也会变为近似值。所以当A和b有微小变化时，即使求解过程精确进行，所得的解相对于原方程组也可能会产生很大的相对误差。对于本节的学习掌握的不是很好，虽然在课后习题中对课堂知识有了一定的巩固，但整体感觉没有很好的掌握它。第四节的迭代法，初次接触迭代法，了解到迭代法就是构造一个无线的向量序列，使他的极限是方程组的解向量。迭代法应考虑收敛性与精度控制的问题。三种迭代方法的基本思想我已经掌握了，但是在matlab 的编程中还存在很大的问题。 在本节的学习中我认为我最大的问题还是程序的编写。通过这段时间的练习，虽然掌握了一些编写方法和技巧。相比于第一章是对其的应用熟练了不少，但在程序编写上还存在很多问题。希望在以后的学习中能尽快熟练掌握它，充分发挥它强大的作用。 二、本章知识梳理 2.1、Gauss消去法（次重点） Gauss消去法基本思想：由消元和回代两个过程组成。 a(k=1,2,```,n-1)均不为零的充分必要条件定理顺序Gauss消去法的前n-1个主元素)(k kk 是方程组的系数矩阵A的前n-1个顺序主子式 题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵 20010-1-24A=0-2131 43 1?? ? ? ? ??? 为例进行具体的求解计算。 一、 算法分析： 一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。 二、 算法设计： 1、 原始实矩阵A R n n ?∈拟上三角化 为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’ 记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于r=1,2,…,n -2执行 (1) 若() (2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。 (2) 计算 令 ()2 ()() 1,1,1,sgn(0,sgn()=1) r r r r r r r r r r r c a a a ρ+++=-=，（若则取 (3) 令-0=r n r u u ?? ? ?? ，()()()-1,2,1(,,...,)r r r T n r r r r r r nr u a c a a ρ++=- (4) 计算 (r)(r)(r)T n-r r+1,r r r+2,r nr r T n-r T n-r n-r n-r n-r r+1r 1u = (a -c ,a ,...,a )ρ I H =I -2uu =H H =I -2u u A =HA H ?? ? ?? (5) 继续 算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。 2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值 记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵， 第一章 题12 给定节点01x =-，11x =，23x =，34x =，试分别对下列函数导出拉格朗日插值余项： （1） （1） 3 ()432f x x x =-+ （2） （2） 4 3 ()2f x x x =- 解 （1）(4) ()0f x =， 由拉格朗日插值余项得(4)0123() ()()()()()()0 4!f f x p x x x x x x x x x ξ-=----=； （2）(4) ()4!f x = 由拉格朗日插值余项得 01234! ()()()()()() 4! f x p x x x x x x x x x -= ----(1)(1)(3)(4)x x x x =+---. 题15 证明：对于()f x 以0x ，1x 为节点的一次插值多项式()p x ，插值误差 012 10()()()max () 8x x x x x f x p x f x ≤≤-''-≤. 证 由拉格朗日插值余项得 01() ()()()()2!f f x p x x x x x ξ''-= --，其中01x x ξ≤≤， 01 0101max ()()()()()()()() 2!2!x x x f x f f x p x x x x x x x x x ξ≤≤''''-=--≤-- 01210()max () 8x x x x x f x ≤≤-''≤. 题22 采用下列方法构造满足条件(0)(0)0p p '==，(1)(1)1p p '==的插值多项式 ()p x ： （1） （1） 用待定系数法； （2） （2） 利用承袭性，先考察插值条件(0)(0)0p p '==，(1)1p =的插值多项式 ()p x . 解 （1）有四个插值条件，故设230123()p x a a x a x a x =+++，2 123()23p x a a x a x '=++， 代入得方程组001231123010231 a a a a a a a a a =? ?+++=?? =? ?++=? 解之，得01230 021 a a a a =??=?? =??=-?数值分析课后题答案
北航数值分析大作业一
数值分析复习题要答案
数值分析作业答案
数值分析第二章复习与思考题
数值分析习题集及答案[1].(优选)
北航数值分析大作业第二题精解
数值分析作业思考题汇总
郑州大学数值分析重点考察内容及各章习题
上海大学_王培康_数值分析大作业
数值分析思考题[综合]
数值分析习题
数值分析第二次大作业
数值分析第二章思考题
北航数值分析大作业第二题
数值分析最佳习题(含答案)
数值分析第二章小结
数值分析大作业QR分解
数值分析课后习题答案
相关主题
文本预览