数值分析作业

数值分析第一章作业1.数值计算方法设计的基本手段是( ).(A) 近似 (B) 插值 (C) 拟合 (D) 迭代2.为了在有限时间内得到结果,用有限过程取代无限过程所产生的近似解与精确解之间的误差称为( ).(A) 舍入误差 (B) 截断误差 (C) 测量误差 (D) 绝对误差3.由于计算机的字长有限,原始数据在机器内的表示以及进行算术运算所产生的误差统称为( ).(A) 舍入误差 (B) 截断误差 (C) 相对误差 (D) 绝对误差4.数值计算方法研究的核心问题可以概括为( )对计算结果的影响.(A) 算法的稳定性 (B) 算法的收敛性 (C) 算法的复杂性 (D) 近似5.当N 充分大时,利用下列各式计算121N N dx I x+=+⎰,等式( )得到的结果最好. (A) arctan(1)arctan()I N N =+- (B) 2arctan(1)I N N =++ (C) 21arctan()1I N N =++ (D) 211I N =+6.计算61), 1.4≈,利用下列哪个公式得到的结果最好?为什么?(B) 3(3- (D) 99-7.计算球体的体积,已知半径的相对误差限不超过3310-⨯,则计算所得体积的相对误差限如何估计?8.设0x >,近似值*x 的相对误差限为δ,试估计*ln x 的误差限.9.计算圆柱体的体积,已知底面半径r 及圆柱高h 的相对误差限不超过δ,则计算所得体积的相对误差限如何估计?.10.用秦九韶算法求32()431f x x x x =-+-在2x =处的值.11.已知近似值 1.0000x *=的误差限4()110x ε*-=⨯,21()16f x x =,求(())f x ε*,并说明x *及()f x *的各有几位有效数字.12.设a 为非零常数,已知0y 的近似值0y *,由递推式1n n y ay -=计算序列{}n y 的近似值,分析该算法的稳定性.。

数值分析大作业一一、算法设计方案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];}}}}三、程序结果。

数值分析上机作业（一）一、算法的设计方案1、幂法求解λ1、λ501幂法主要用于计算矩阵的按模最大的特征值和相应的特征向量，即对于|λ1|≥|λ2|≥.....≥|λn|可以采用幂法直接求出λ1，但在本题中λ1≤λ2≤……≤λ501，我们无法判断按模最大的特征值。

但是由矩阵A的特征值条件可知|λ1|和|λ501|之间必然有一个是最大的，通过对矩阵A使用幂法迭代一定次数后得到满足精度ε=10−12的特征值λ0，然后在对矩阵A做如下的平移：B=A-λ0I由线性代数(A-PI)x=(λ-p)x可得矩阵B的特征值为：λ1-λ0、λ2-λ0…….λ501-λ0。

对B矩阵采用幂法求出B矩阵按模最大的特征值为λ∗=λ501-λ0，所以λ501=λ∗+λ0，比较λ0与λ501的大小，若λ0>λ501则λ1=λ501，λ501=λ0；若λ0<λ501，则令t=λ501，λ1=λ0，λ501=t。

求矩阵M按模最大的特征值λ的具体算法如下：任取非零向量u0∈R nηk−1=u T(k−1)∗u k−1y k−1=u k−1ηk−1u k=Ay k−1βk=y Tk−1u k(k=1,2,3……)当|βk−βk−1||βk|≤ε=10−12时，迭终终止，并且令λ1=βk2、反幂法计算λs和λik由已知条件可知λs是矩阵A 按模最小的特征值，可以应用反幂法直接求解出λs。

使用带偏移量的反幂法求解λik，其中偏移量为μk=λ1+kλ501−λ140(k=1,2,3…39)，构造矩阵C=A-μk I，矩阵C的特征值为λik−μk，对矩阵C使用反幂法求得按模最小特征值λ0，则有λik=1λ0+μk。

求解矩阵M按模最小特征值的具体算法如下：任取非零向量u 0∈R n ηk−1= u T (k−1)∗u k−1y k−1=u k−1ηk−1 Au k =y k−1βk =y T k−1u k (k=1,2,3……)在反幂法中每一次迭代都要求解线性方程组Au k =y k−1，当K 足够大时，取λn =1βk 。

问题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代表端点的一阶导。

数值分析大作业（二）院系：学号：姓名：一、算法设计方案：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. 误差分为有__模型误差___, _观测误差___, __方法误差____, ___舍入误差____, Taylor 展开式近似表达函数产生的误差是_方法误差 .2. 插值余项是插值多项式的 方法误差。

0.2499作为1/4的近似值, 有几位有效数字？00.24990.249910,0m =⨯=即，031|0.2499|0.00010.5100.510,34m n n ---=<⨯=⨯=即22 3.1428751...,7=作为圆周率的近似值，误差和误差限分别是多少，有几位有效数字？2133.142875 3.14159260.00126450.5100.510---=<⨯=⨯有3位有效数字.* 有效数字与相对误差的关系3. 利用递推公式计算积分110,1,2,...,9n x n I x e dx n -==⎰错误!未找到引用源。

, 建立稳定的数值算法。

该算法是不稳定的。

因为:11()()...(1)!()n n n I n I n I εεε-=-==-111n n I I n n -=-， 10110I =4. 衡量算法优劣的指标有__时间复杂度,__空间复杂度_.时间复杂度是指: , 两个n 阶矩阵相乘的乘法次数是 , 则称两个n 阶矩阵相乘这一问题的时间复杂度为 .二 代数插值1.根据下表数据建立不超过二次的Lagrange 和Newton 插值多项式, 并写出误差估计式, 以及验证插值多项式的唯一性。

x 0 1 4f(x) 1 9 3Lagrange:设0120120,1,4;()1()9()3x x x f x f x f x ======则，， 对应 的标准基函数 为:1200102()()(1)(x 4)1()(1)(x 4)()()(01)(04)4x x x x x l x x x x x x ----===------ 1()...l x =2()...l x =因此, 所求插值多项式为:220()()()....i i i P x f x l x ===∑ (3)2()()(0)(1)(x 4)3!f R x x x ξ=--- Newton:构造出插商表:xi f(xi ) 一 二 三0 11 9 84 3 -2 -5/2所以, 所求插值多项式为:2001001201()()[,]()[,,]()()518(0)(0)(1)2...P x f x f x x x x f x x x x x x x x x x =+-+--=+----=插值余项: 2()[0,1,4,](0)(1)(x 4)R x f x x x =---2. 已知函数f(0)=1,f(1)=3,f(2)=7,则f[0,1]=___2________, f[0,1,2]=____1______)('],[000x f x x f =3.过0,1两节点构造三次Hermite 插值多项式, 使得满足插值条件: f(0)=1. .’(0)=... f(1.=2. .’(1)=1设0101010,1,()1()2'()0,'()1x x f x f x f x f x ======则，， 写出插商表:xi f(xi) 一 二 三0 10 1 01 a 1 11 a 1 0 a-1因此, 所求插值多项式为:插值余项:222()[0,0,1,1,](1)R x f x x x =-4.求f(x)=sinx 在[a,b]区间上的分段线性插值多项式, 并写出误差估计式。

数值分析大作业数值分析大作业姓名：黄晨晨学号：S1*******学院：储运与建筑工程学院学院班级：储建研17-2实验3.1 Gauss消去法的数值稳定性实验实验目的：理解高斯消元过程中出现小主元即很小时引起方程组解数值不定性实验内容：求解方程组Ax=b，其中（1）A1=0.3×10?1559.14315.291?6.130?1211.29521211，b1=59.1746.7812;（2）A2=10?7013 2.099999999999625?15?10102，b2=85.90000000000151;实验要求：(1)计算矩阵的条件数，判断系数矩阵是良态的还是病态的(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4(4)观察小主元并分析对计算结果的影响（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的代码：format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2]n=4C1=cond(A1,1) %C1为A1矩阵1范数下的条件数C2=cond(A1,2) %C2为A1矩阵2范数下的条件数C3=cond(A1,inf) %C3为1矩阵谱范数下的条件数结果：C1 =1.362944708720448e+02C2 =6.842955771253409e+01C3 =8.431146*********e+01显然A1矩阵为病态矩阵将矩阵A2，b2输入上述代码中求得A2矩阵的条件数为：C1 =1.928316831682894e+01C2 =8.993938090170119e+00C3 =1.835643564356072e+01显然A2矩阵也为病态矩阵(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4代码：for k=1:n-1a=max(abs(A1(k:n,k)))if a==0returnendfor i=k:nif abs(A1(i,k))==ay=A1(i,:)A1(i,:)=A1(k,:)A1(k,:)=yx=b1(i,:)b1(i,:)=b1(k,:)b1(k,:)=xbreakendendif A1(k,k)~=0A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k)A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n) elsebreakendendL=tril(A1,0);for i=1:nL(i,i)=1;endLU=triu(A1,0)y1=L\b1x1=U\y1得到如下结果：x1 =3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2]b2=[8;5.900000000001;5;1]代入上述代码求得结果如下：X2 =4.440892098500626e-16-9.999999999999993e-019.999999999999997e-011.000000000000000e+00(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4代码：format longeformat compactA1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1] b1=[59.17;46.78;1;2][L,U]=lu(A1)y1=L\b1x1=U\y1求得如下结果：x1=3.845714853511634e+001.609517394778522e+00-1.547605454206655e+011.041130489899787e+01将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2] b2=[8;5.900000000001;5;1]代入上述代码，求得结果如下：x 2 =4.440892098500626e-16 -9.999999999999993e-01 9.999999999999997e-01 9.999999999999999e-01（2）（3）求得结果相同，可验证结果正确。

数值分析课后作业：习题一1.在字长为3的十进制计算机上计算f （3.33）和g （3.33），其中f(x)=x 4-x 3+3x 2+x-2,g(x)=(((x-1)x+3)x+1)x-2解： m=3; f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m); g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m); f(3.33) g(3.33) 有ans = 121 ans =121 2.下列各近似值的绝对误差限都是1021⨯-3，试指出它们各有几位有效数字：x=1.00052， y=0.05， z=0.00052.解：当 x=1.00052时， 由丨X*—X 丨 ≤0.5×10-3 得 x=1.00052 有四位有效数字； 同理 y=0052 有两位有效数字 Z=0.00052有零位有效数字 3,计算圆的面积，要使其相对误差限为1%，问测量半径r 允许的相对误差限是多少？ 解：设圆的面积为S ， 由题意有|e(S)|≤1%。

又S=πr 2 dS=2πr dr 所以 dS/S=(2πrdr)/(πr 2)=2(dr/r)∴|e(r)|≈21|e(S)|≤0.5×1%=0.5% 11.数组与矩阵是Matlab 编程的基础，试学习Matlab 的数组与矩阵的表示方法，并举例介绍数组、矩阵的常见运算. 解：>> syms a b c d; >> a=[1 2 3];>> b=[4 5 6];>> a+bans =5 7 9>> b-aans =3 3 3>> a.*bans =4 10 18 >> a.^2 ans = 1 4 9>> c=[1 2 3;1 2 3;1 2 3];>> d=[4 5 6;4 5 6;4 5 6];>> cc = 1 2 3 1 2 3 1 2 3d = 4 5 6 4 5 6 4 5 6 >> c+dans =5 7 9 5 7 9 5 7 9>> d-cans = 3 3 33 3 33 3 3 12.学习使用Matlab 命令help 和doc 学习自己感兴趣的Matlab 的运算、函数或命令的用法，并对于任意给定的实数a,b,c,编写Matlab 程序求方程ax 2+bx+c=0的根. 解：x 1=a ac b b b 24)sgn(2---, x 2=1ax c1 x>0 其中 sgn = 0 x=0 -1 x<0 disp('Please input the coefficients of');disp('quadratic equation ax^2+bx+c=0, respectively') a=input('a='); b=input('b='); c=input('c=');m=3; if abs(a)<eps & abs(b)<eps error End if abs(a)<eps disp('Since a=0, quadrtic equation degen erates into a linear equation.') disp('The only solution of the linear equtio n is')x=digit(-c/b,m) return Enddelta=b^2-4*a*c; temp=sqrt(delta); x 1=(-b+temp)/(2*a) ; x 2=(-b-temp)/(2*a) ;err1=abs(a*x 1^2+b*x 1+c) ; err2=abs(a*x 2^2+b*x 2+c) ; if b>0x 1=(-b-temp)/(2*a) End if b<0x 1=(-b+temp)/(2*a) End if b=0x 1=temp/(2*a) Endx 2=c/(a*x 1)err1=abs(a*x 1^2+b*x 1+c) err2=abs(a*x 2^2+b*x 2+c) if abs(a)<epsdisp('Since a=0, quadrtic equation degen erates into a linear equation.')disp('The only solution of the linear equtio n is')x=digit(-c/b,m) return Enddelta=digit(digit(b^2,m)-digit(4*digit(a*c,m),m),m);temp=digit(sqrt(delta),m);x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); x 2=digit(digit(-b-temp,m)/digit(2*a,m),m); err1=abs(a*x 1^2+b*x 1+c); err2=abs(a*x 2^2+b*x 2+c); if b>0x 1=digit(digit(-b-temp,m)/digit(2*a,m),m) ; End if b<0x 1=digit(digit(-b+temp,m)/digit(2*a,m),m); End if b=0x 1=digit(temp/digit(2*a,m),m); Endx 2=digit(digit(c/a,m)/x1,m) ; err1=abs(a*x 1^2+b*x 1+c) ; err2=abs(a*x 2^2+b*x 2+c) ; 14分别利用ln (1+x)=11,)1(11≤<--+∞=∑x nx nn n 和ln11...),12...53(2111253<<-++++++=-++x n x x x x x x n ，给出计算ln2的近似方法，编写相应的Matlab 程序，并比较算法运行情况. 解：方法一： x=1; s=0;for k=1:100s=s+(-1)^(k+1)*(x^k)/k; end sq=log(2)err=abs(t-q) ans= t =0.6882 q =0.6931 err = 0.0050方法二x=1/3; s=0;for k=1:2:100 s=s+(x^k)/k; end t=2*s q=log(2)err=abs(t-q) Ans= t =0.6931 q =0.6931 err =2.2204e-16所以方法二较方法一好。