数值分析大作业模板

数值分析大作业三四五六七数值分析大作业三四五六七Document number【SA80SAB-SAA9SYT-SAATC-SA6UT-SA18】大作业三1. 给定初值0x 及容许误差，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下：函数m 文件：fu.mfunction Fu=fu(x)Fu=x^3/3-x;end函数m 文件：dfu.mfunction Fu=dfu(x)Fu=x^2-1;end用Newton 法求根的通用程序Newton.mclear;x0=input('请输入初值x0:');ep=input('请输入容许误差:');flag=1;while flag==1x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<ep< p="">flag=0;endx0=x1;endfprintf('方程的一个近似解为：%f\n',x0);寻找最大δ值的程序：Find.mcleareps=input('请输入搜索精度：');ep=input('请输入容许误差：');flag=1;k=0;x0=0;while flag==1sigma=k*eps;x0=sigma;k=k+1;m=0;flag1=1;while flag1==1 && m<=10^3x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)endm=m+1;x0=x1;endif flag1==1||abs(x0)>=epflag=0;endendfprintf('最大的sigma 值为：%f\n',sigma);2．求下列方程的非零根5130.6651()ln 05130.665114000.0918x x f x x +??=-= ?-解：Matlab 程序为：（1）主程序clearclcformat longx0=765;N=100;errorlim=10^(-5);x=x0-f(x0)/subs(df(),x0);n=1;while n<n< p="">x=x0-f(x0)/subs(df(),x0);if abs(x-x0)>errorlimn=n+1;elsebreak;endx0=x;enddisp(['迭代次数: n=',num2str(n)])disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)])（2）子函数非线性函数ffunction y=f(x)y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918);end（3）子函数非线性函数的一阶导数dffunction y=df()syms x1y=log((513+0.6651*x1)/(513-0.6651*x1))-x1/(1400*0.0918);y=diff(y);end运行结果如下：迭代次数: n=5所求非零根: 正根x1=767.3861 负根x2=-767.3861大作业四试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式. 对函数21()14f x x=+在区间[-5，5]上实现10次多项式插值.分析：（1）输出插值多项式。

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

一、题目：关于x, y, t, u, v, w 的下列方程组0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩1、试用数值方法求出f(x, y)在区域 {(,)|00.8,0.5 1.5}D x y x y =≤≤≤≤上的一个近似表达式,0(,)kr s rsr s p x y cx y ==∑要求(,)p x y 一最小的k 值达到以下的精度10202700((,)(,))10i j i j i j f x y p x y σ-===-≤∑∑其中，0.08,0.50.05i j x i y j ==+。

2、计算****(,),(,)i j i j f x y p x y （i = 1, 2, …，8；j = 1, 2,…，5）的值，以观察(,)p x y 逼近(,)f x y 的效果，其中，*i x =0.1i , *j y =0.5+0.2j 。

说明：1、用迭代方法求解非线性方程组时，要求近似解向量()k x 满足()(1)()12||||/||||10k k k x x x --∞∞-≤2、作二元插值时，要使用分片二次代数插值。

3、要由程序自动确定最小的k 值。

4、打印以下内容：●算法的设计方案。

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

●数表：,,i j x y (,)i j f x y （i = 0,1,2,…,10；j = 0,1,2,…,20）。

●选择过程的,k σ值。

●达到精度要求时的,k σ值以及(,)p x y 中的系数rs c （r = 0,1,…,k;s = 0,1,…,k ）。

●数表：**,,i j x y ****(,),(,)i j i j f x y p x y （i = 1, 2, ...,8；j = 1, 2, (5)。

数值分析上机作业（一）一、算法的设计方案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 。

《数值分析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,j2.求解λ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<stdio.h>#include<iostream.h>#include<stdlib.h>#include<math.h>#include<float.h>#include<iomanip.h>#include<time.h>#define E 1.0e-12 /*定义全局变量相对误差限*/int max2(int a,int b) /*求两个整型数最大值的子程序*/{if(a>b)return a;elsereturn b;}int min2(int a,int b) /*求两个整型数最小值的子程序*/{if(a>b)return b;elsereturn a;}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);}void assignment(double array[5][501]) /*将矩阵A转存为数组C[5][501]*/{int i,j,k;//所有元素归零for(i=0;i<=4;){for(j=0;j<=500;){array[i][j]=0;j++;}i++;}//第0,4行赋值for(j=2;j<=500;){k=500-j;array[0][j]=-0.064;array[4][k]=-0.064;j++;}//第1,3行赋值for(j=1;j<=500;){k=500-j;array[1][j]=0.16;array[3][k]=0.16;j++;}//第2行赋值for(j=0;j<=500;){ k=j;j++;array[2][k]=(1.64-0.024*j)*sin((double)(0.2*j))-0.64*exp((double)(0.1/j));}}double mifa(double u[501],double array[5][501],double p) /*带原点平移的幂法*/ {int i,j; /* u[501]为初始迭代向量*/double a,b,c=0; /* array[5][501]为矩阵A的转存矩阵*/double y[501]; /*p为平移量*/for(;;){a=0;b=0;/*选用第一种迭代格式*///求ηk-1for(i=0;i<=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求y k-1for(i=0;i<=500;i++){y[i]=u[i]/a;}//求u kfor(i=0;i<=500;i++){u[i]=0;for(j=max2(i-2,0);j<=min2(i+2,500);j++){u[i]+=array[i-j+2][j]*y[j];}u[i]=u[i]-p*y[i]; /*引入平移量*/}//求βkfor(i=0;i<=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)<=E) /*达到精度水平，迭代终止*/break;c=b;}return (b+p); /*直接返回A的特征值*/}void chuzhi(double a[]) /*用随机数为初始迭代向量赋值*/ {int i;srand((int)time(0));for(i=0;i<=500;i++){a[i]=(10.0*rand()/RAND_MAX); /*生成0~10的随机数*/}}void chuzhi2(double a[],int j) /*令初始迭代向量为e i*/{int i;for(i=0;i<=500;i++){a[i]=0;}a[j]=1;}void LU(double array[5][501]) /*对矩阵A进行Doolittle分解*/{ /*矩阵A转存在C[5][501]中*/int j,k,t; /*分解结果L，U分别存在C[5][501]的上半部与下半部*/ for(k=0;k<=500;k++){for(j=k;j<=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t<=(k-1);t++){array[k-j+2][j]-=array[k-t+2][t]*array[t-j+2][j];}}if(k<500)for(j=k+1;j<=min2((k+2),500);j++){for(t=max3(0,k-2,j-2);t<=(k-1);t++){array[j-k+2][k]-=array[j-t+2][t]*array[t-k+2][k];}array[j-k+2][k]=array[j-k+2][k]/array[2][k];}}}double fmifa(double u[501],double array[5][501],double p){ /*带原点平移的反幂法*/ int i,j;double a,b,c=0;double y[501];//引入平移量for(i=0;i<=500;i++){array[2][i]-=p;}//先将矩阵Doolittle分解LU(array);for(;;){a=0;b=0;//求ηk-1for(i=0;i<=500;i++){a=a+u[i]*u[i];}a=sqrt(a);//求y k-1for(i=0;i<=500;i++){y[i]=u[i]/a;}//回带过程，求解u kfor(i=0;i<=500;i++){u[i]=y[i];}for(i=1;i<=500;i++){for(j=max2(0,(i-2));j<=(i-1);j++){u[i]-=array[i-j+2][j]*u[j];}}u[500]=u[500]/array[2][500];for(i=499;i>=0;i--){for(j=i+1;j<=min2((i+2),500);j++){u[i]-=array[i-j+2][j]*u[j];}u[i]=u[i]/array[2][i];}//求βkfor(i=0;i<=500;i++){b+=y[i]*u[i];}if(fabs((b-c)/b)<=E) /*达到精度要求，迭代终止*/break;c=b;}return (p+(1/b)); /*直接返回距离原点P最接近的A的特征值*/ }//主函数main(){ int i;double d1,d501,ds,d,a;double u[501];double MatrixC[5][501];printf(" 《数值分析》计算实习题目第一题\n");printf(" SY1103120 朱舜杰\n");//将矩阵A转存为MatrixCassignment(MatrixC);//用带原点平移的幂法求解λ1，λ501chuzhi(u);d=mifa(u,MatrixC,0);chuzhi(u);a=mifa(u,MatrixC,d);if(d<0){d1=d;d501=a;}else{d501=d;d1=a;}printf("λ1=%.12e\n",d1);printf("λ501=%.12e\n",d501);//用反幂法求λschuzhi(u);ds=fmifa(u,MatrixC,0);printf("λs=%.12e\n",ds);//用带原点平移的反幂法求λikfor(i=1;i<=39;i++){a=d1+(i*(d501-d1))/40;assignment(MatrixC);chuzhi(u);d=fmifa(u,MatrixC,a);printf("与μ%02d=%+.12e最接近的特征值λi%02d=%+.12e\n",i,a,i,d);}//求A的条件数d=fabs((d1/ds));printf("A的（谱范数）条件数cond<A>2=%.12e\n",d);//求detAassignment(MatrixC);LU(MatrixC);a=1;for(i=0;i<=500;i++){a*=MatrixC[2][i];}printf("行列式detA=%.12e\n",a);//测试不同迭代初始向量对λ1计算结果的影响。

数值分析大作业（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的整数位数） 所以有绝对误差限 14311(1)101022x ε--≤⨯=⨯相对误差限 31()0.510(1)0.00923%5.4201r x x x εε-⨯=== 对x2：有a2=5, m2=0 所以有绝对误差限 04411(2)101022x ε--≤⨯=⨯相对误差限 42()0.510(2)0.00923%0.54202r x x x εε-⨯=== 对x3：有a3=5, m3=-2 所以有绝对误差限 23511(3)101022x ε---≤⨯=⨯相对误差限 53()0.510(3)0.0923%0.005423r x x x εε-⨯=== 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022x ε-≤⨯=相对误差限 4()0.5(4)0.0083%60004r x x x εε=== 对x5：有a5=6, m5=5 所以有绝对误差限 51411(5)101022x ε-≤⨯=⨯相对误差限 45()0.510(5)8.3%600005r x x x εε⨯===2.对矩阵A 进行LU 分解， 并求解方程组Ax b =其中211132122A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦，465b ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦解：A=LU 代入方程Ax b = 可转化为L y bU x y⎧=⎪⎨=⎪⎩先对矩阵A 进行LU 分解，如下11121311121321222321112112222113233132333111311232223113322333111u u u u u u A LU l u u l u l u u l u u l l u l u l u l u l u l u u ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥===++⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥+++⎣⎦⎣⎦⎣⎦根据系数相应相等有：第一行：112u =，121u =，131u = 第二行：21112121l u l ==，可得210.5l =211222220.513l u u u +=⨯+=，可得22 2.5u = 211323230.512l u u u +=⨯+=，可得23 1.5u =第三行：31113121l u l ==，可得310.5l =31123222320.51 2.52l u l u l +=⨯+=，可得320.6l =3113322333330.510.6 1.52l u l u u u ++=⨯+⨯+=，可得330.6u =所以有：12110.51 2.5 1.50.50.610.6A LU ⎡⎤⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦解方程如下123140.5160.50.615y y y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦，可得123440.6y y y ⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦ 11223321142.5 1.540.60.6x y x y x y ⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦，可得123111x x x ⎡⎤⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦3. 用 J 迭代法和 G-S 迭代法求解方程组 123123123202324812231530x x x x x x x x x ++=⎧⎪++=⎨⎪-+=⎩时， 若取初始解向量(0)(0,0,0)T x = ， 问各需迭代多少次才能使误差()*610k x x-∞-≤ 。

北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号ZY*******学生姓名许阳教师孙玉泉日期2021 年11月26 日设有501501⨯的实对称矩阵A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=5011A a b c b c c b c b a其中，064.0,16.0),501,,2,1(64.0)2.0sin()024.064.1(1.0-==⋅⋅⋅=--=c b i e i i a ii 。

矩阵A 的特征值为)501,,2,1(⋅⋅⋅=i i λ，并且有||min ||,501150121i i s λλλλλ≤≤=≤⋅⋅⋅≤≤1λ，501λ和s λ的值。

A 的与数4015011λλλμ-+=kk 最接近的特征值)39,,2,1(⋅⋅⋅=k k i λ。

A 的(谱范数)条件数2)A (cond 和行列式detA 。

一 方案设计1 求1λ，501λ和s λ的值。

s λ为按模最小特征值，||min ||5011i i s λλ≤≤=。

可使用反幂法求得。

1λ，501λ分别为最大特征值及最小特征值。

可使用幂法求出按模最大特征值，如结果为正，即为501λ，结果为负，那么为1λ。

使用位移的方式求得另一特征值即可。

2 求A 的与数4015011λλλμ-+=kk 最接近的特征值)39,...,2,1(=k k i λ。

题目可看成求以k μ为偏移量后，按模最小的特征值。

即以k μ为偏移量做位移，使用反幂法求出按模最小特征值后，加上k μ，即为所求。

3 求A 的(谱范数)条件数2)(A cond 和行列式detA 。

矩阵A 为非奇异对称矩阵，可知，||)(min max2λλ=A cond(1-1)其中m ax λ为按模最大特征值，min λ为按模最小特征值。

detA 可由LU 分解得到。

因LU 均为三角阵，那么其主对角线乘积即为A 的行列式。

二 算法实现1 幂法使用如下迭代格式：⎪⎪⎩⎪⎪⎨⎧⋅===⋅⋅⋅=------||max |)|sgn(max ||max /),,(111111)0()0(10k k k k k k k k Tn u u Ay u u u y u u u β任取非零向量 (2-1)终止迭代的控制理论使用εβββ≤--||/||1k k k ， 实际使用εβββ≤--||/||||||1k k k(2-2)由于不保存A 矩阵中的零元素，只保存主对角元素a[501]及b,c 值。

《数值分析》大作业四一、算法设计方案：复化梯形积分法，选取步长为1/500=0.002，迭代误差控制在E ≤1.0e-10①复化梯形积分法：11()[()()2()]2n bak hf x dx f a f b f a kh -=⎰≈+++∑，截断误差为：322()''()''(),[,]1212T b a b a R f h f a b n ηηη--=-=-∈其中。

复化Simpson 积分法，选取步长为1/50=0.02，迭代误差控制在E ≤1.0e-10②Simpson 积分法：121211()[()()4()2()]3m m bi i a i i hf x dx f a f b f x f x --==≈+++∑∑⎰， 截断误差为：4(4)(),[,]180s b a R h f a b ηη-=-∈。

③Guass积分法选用Gauss-Legendre 求积公式：111()()ni i i f x dx A f x -=≈∑⎰截断误差为：R= ()()n 2n 422n!2×(2[2!]2n 1f n n ⨯（2）η（））＋ η∈(1,1)。

选择9个节点：-0.9681602395,-0.8360311073,-0.6133714327,-0.3242534234,0,0.3242534234,0.6133714327,0.8360311073,0.9681602395， 对应的求积系数依次为：0.0812743884,0.1806481607,0.2606106964,0.3123470770,0.3302393550,0.3123470770,0.2606106964,0.1806481607,0.0812743884。

二、程序源代码：#include<stdio.h>#include<math.h>#include<stdlib.h>#define E 1.0e-10/****定义函数g和K*****/double g(double a){double b;b=exp(4*a)+(exp(a+4)-exp(-a-4))/(a+4);return b;}double K(double a,double b){double c;c=exp(a*b);return c;}/******复化梯形法******/void Tixing( ){double u[1001],x[1001],h,c[1001],e;int i,j,k;FILE *fp;fp=fopen("f:/result0. xls ","w");h=1.0/1500;for(i=0;i<3001;i++){x[i]=i*h-1;u[i]=g(x[i]);}for(k=0;k<100;k++){e=0;for(i=0;i<1001;i++){for(j=1,c[i]=0;j<N-1;j++)c[i]+=K(x[i],x[j])*u[j];u[i]=g(x[i])-h*c[i]-h/2*(K(x[i],x[0])*u[0]+K(x[i],x[N-1])*u[N-1]);e+=h*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);}if(e<=E) break;}for(i=0;i<1001;i++)fprintf(fp,"%.12lf,%.12lf\n",x[i],u[i]);fclose(fp);}/******复化Simpson法******/void simpson( ){double u[101],x[101],h,c[101],d[101],e;int i,j,k;FILE *fp;fp=fopen("f:/result1.xls","w");h=1.0/50;for(i=0;i<101;i++){x[i]=i*h-1;u[i]=g(x[i]);}for(k=0;k<50;k++){e=0;for(i=0;i<101;i++){for(j=1,c[i]=0,d[i]=0;j<51;j++){c[i]+=K(x[i],x[2*j-1])*u[2*j-1];if(j<50)d[i]+=K(x[i],x[2*j])*u[2*j];}u[i]=g(x[i])-4*h/3*c[i]-2*h/3*d[i]-h/3*(K(x[i],x[0])*u[0]+K(x[i],x[M-1])*u[M-1]);e+=h*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);}if(e<=E) break;}for(i=0;i<101;i++)fprintf(fp,"%.12lf,%.12lf\n",x[i],u[i]);fclose(fp);}/******Gauss积分法******/void gauss( ){double x[9]={-0.9681602395,-0.8360311073,-0.6133714327,-0.3242534234,0,\0.3242534234,0.6133714327,0.8360311073,0.9681602395},A[9]={0.0812743884,0.1806481607,0.2606106964,0.3123470770,0.3302393550,\0.3123470770,0.2606106964,0.1806481607,0.0812743884},u[9],c[9],e;int i,j,k;FILE *fp;fp=fopen("f:/result2. xls ","w");for(i=0;i<9;i++)u[i]=g(x[i]);for(k=0;k<50;k++){e=0;for(i=0;i<9;i++){for(j=0,c[i]=0;j<9;j++)c[i]+=A[j]*K(x[i],x[j])*u[j];u[i]=g(x[i])-c[i];e+=A[i]*(exp(4*x[i])-u[i])*(exp(4*x[i])-u[i]);}if(e<=E) break;}for(i=0;i<9;i++)fprintf(fp,"%.12lf,%.12lf\n",x[i],u[i]);fclose(fp);}/******主函数******/main(){Tixing ( );Simpson( );Gauss( );return 0;}三、运算结果复化梯形数据-10.018323-0.920.02523-0.9980.018471-0.9180.025433-0.9960.018619-0.9160.025637-0.9940.018768-0.9140.025843-0.9920.018919-0.9120.026051-0.990.019071-0.910.02626-0.9880.019224-0.9080.026471-0.9860.019378-0.9060.026683-0.9840.019534-0.9040.026897-0.9820.019691-0.9020.027113-0.980.019849-0.90.027331-0.9780.020008-0.8980.02755-0.9760.020169-0.8960.027772-0.9740.020331-0.8940.027995-0.9720.020494-0.8920.028219-0.970.020658-0.890.028446-0.9680.020824-0.8880.028674-0.9660.020992-0.8860.028905-0.9640.02116-0.8840.029137-0.9620.02133-0.8820.029371-0.960.021501-0.880.029607-0.9580.021674-0.8780.029844-0.9560.021848-0.8760.030084-0.9540.022023-0.8740.030326-0.9520.0222-0.8720.030569-0.950.022378-0.870.030815-0.9480.022558-0.8680.031062-0.9460.022739-0.8660.031311-0.9440.022922-0.8640.031563-0.9420.023106-0.8620.031816-0.940.023291-0.860.032072-0.9380.023478-0.8580.032329-0.9360.023667-0.8560.032589-0.9340.023857-0.8540.032851-0.9320.024048-0.8520.033114-0.930.024241-0.850.03338-0.9280.024436-0.8480.033648-0.9260.024632-0.8460.033918-0.9240.02483-0.8440.034191-0.9220.025029-0.8420.034465-0.840.034742-0.760.047841-0.8380.035021-0.7580.048225-0.8360.035302-0.7560.048613 -0.8340.035586-0.7540.049003 -0.8320.035872-0.7520.049396 -0.830.03616-0.750.049793 -0.8280.03645-0.7480.050193 -0.8260.036743-0.7460.050596 -0.8240.037038-0.7440.051002 -0.8220.037335-0.7420.051412 -0.820.037635-0.740.051825 -0.8180.037937-0.7380.052241 -0.8160.038242-0.7360.052661 -0.8140.038549-0.7340.053084 -0.8120.038858-0.7320.05351 -0.810.039171-0.730.05394 -0.8080.039485-0.7280.054373 -0.8060.039802-0.7260.054809 -0.8040.040122-0.7240.05525 -0.8020.040444-0.7220.055693 -0.80.040769-0.720.056141 -0.7980.041096-0.7180.056591 -0.7960.041426-0.7160.057046 -0.7940.041759-0.7140.057504 -0.7920.042094-0.7120.057966 -0.790.042432-0.710.058431 -0.7880.042773-0.7080.058901 -0.7860.043116-0.7060.059374 -0.7840.043463-0.7040.05985 -0.7820.043812-0.7020.060331 -0.780.044164-0.70.060816 -0.7780.044518-0.6980.061304 -0.7760.044876-0.6960.061796 -0.7740.045236-0.6940.062293 -0.7720.045599-0.6920.062793 -0.770.045966-0.690.063297 -0.7680.046335-0.6880.063805 -0.7660.046707-0.6860.064318 -0.7640.047082-0.6840.064834 -0.7620.04746-0.6820.065355-0.680.06588-0.60.090722 -0.6780.066409-0.5980.091451-0.6760.066942-0.5960.092185 -0.6740.06748-0.5940.092926 -0.6720.068022-0.5920.093672 -0.670.068568-0.590.094424 -0.6680.069119-0.5880.095183 -0.6660.069674-0.5860.095947 -0.6640.070234-0.5840.096718 -0.6620.070798-0.5820.097494 -0.660.071366-0.580.098277 -0.6580.071939-0.5780.099067 -0.6560.072517-0.5760.099862 -0.6540.0731-0.5740.100664 -0.6520.073687-0.5720.101473 -0.650.074278-0.570.102288 -0.6480.074875-0.5680.103109 -0.6460.075476-0.5660.103937 -0.6440.076082-0.5640.104772 -0.6420.076694-0.5620.105614 -0.640.077309-0.560.106462 -0.6380.07793-0.5580.107317 -0.6360.078556-0.5560.108179 -0.6340.079187-0.5540.109048 -0.6320.079823-0.5520.109924 -0.630.080464-0.550.110806 -0.6280.08111-0.5480.111696 -0.6260.081762-0.5460.112593 -0.6240.082418-0.5440.113498 -0.6220.08308-0.5420.114409 -0.620.083748-0.540.115328 -0.6180.08442-0.5380.116254 -0.6160.085098-0.5360.117188 -0.6140.085782-0.5340.118129 -0.6120.086471-0.5320.119078 -0.610.087165-0.530.120035 -0.6080.087865-0.5280.120999 -0.6060.088571-0.5260.12197 -0.6040.089282-0.5240.12295 -0.6020.089999-0.5220.123938-0.550.110806-0.470.152592 -0.5480.111696-0.4680.153817-0.5460.112593-0.4660.155053-0.5440.113498-0.4640.156298-0.5420.114409-0.4620.157553-0.540.115328-0.460.158819-0.5380.116254-0.4580.160095-0.5360.117188-0.4560.16138-0.5340.118129-0.4540.162677-0.5320.119078-0.4520.163983-0.530.120035-0.450.1653-0.5280.120999-0.4480.166628-0.5260.12197-0.4460.167966-0.5240.12295-0.4440.169315-0.5220.123938-0.4420.170675-0.520.124933-0.440.172046-0.5180.125936-0.4380.173428-0.5160.126948-0.4360.174821-0.5140.127967-0.4340.176225-0.5120.128995-0.4320.17764-0.510.130031-0.430.179067-0.5080.131076-0.4280.180505-0.5060.132128-0.4260.181955-0.5040.13319-0.4240.183416-0.5020.134259-0.4220.18489-0.50.135338-0.420.186375-0.4980.136425-0.4180.187871-0.4960.13752-0.4160.18938-0.4940.138625-0.4140.190901-0.4920.139738-0.4120.192435-0.490.140861-0.410.19398-0.4880.141992-0.4080.195538-0.4860.143132-0.4060.197109-0.4840.144282-0.4040.198692-0.4820.145441-0.4020.200288-0.480.146609-0.40.201897-0.4780.147786-0.3980.203518-0.4760.148973-0.3960.205153-0.4740.15017-0.3940.206801-0.4720.151376-0.3920.208462-0.390.210136-0.310.289382 -0.3880.211824-0.3080.291706-0.3860.213525-0.3060.294049 -0.3840.21524-0.3040.296411 -0.3820.216969-0.3020.298792 -0.380.218711-0.30.301192 -0.3780.220468-0.2980.303611 -0.3760.222239-0.2960.306049 -0.3740.224024-0.2940.308508 -0.3720.225823-0.2920.310985 -0.370.227637-0.290.313483 -0.3680.229465-0.2880.316001 -0.3660.231308-0.2860.318539 -0.3640.233166-0.2840.321098 -0.3620.235039-0.2820.323677 -0.360.236927-0.280.326277 -0.3580.23883-0.2780.328897 -0.3560.240748-0.2760.331539 -0.3540.242682-0.2740.334202 -0.3520.244631-0.2720.336886 -0.350.246596-0.270.339592 -0.3480.248576-0.2680.34232 -0.3460.250573-0.2660.345069 -0.3440.252586-0.2640.347841 -0.3420.254614-0.2620.350635 -0.340.256659-0.260.353451 -0.3380.258721-0.2580.35629 -0.3360.260799-0.2560.359151 -0.3340.262894-0.2540.362036 -0.3320.265005-0.2520.364944 -0.330.267134-0.250.367875 -0.3280.269279-0.2480.37083 -0.3260.271442-0.2460.373809 -0.3240.273622-0.2440.376811 -0.3220.27582-0.2420.379838 -0.320.278035-0.240.382888 -0.3180.280268-0.2380.385964 -0.3160.28252-0.2360.389064 -0.3140.284789-0.2340.392189 -0.3120.287076-0.2320.395339-0.230.398514-0.150.548804-0.2280.401715-0.1480.553212-0.2260.404942-0.1460.557655 -0.2240.408194-0.1440.562134 -0.2220.411473-0.1420.56665 -0.220.414778-0.140.571201 -0.2180.418109-0.1380.575789 -0.2160.421467-0.1360.580414 -0.2140.424853-0.1340.585076 -0.2120.428265-0.1320.589775 -0.210.431705-0.130.594512 -0.2080.435172-0.1280.599287 -0.2060.438668-0.1260.604101 -0.2040.442191-0.1240.608953 -0.2020.445743-0.1220.613844 -0.20.449323-0.120.618774 -0.1980.452932-0.1180.623744 -0.1960.45657-0.1160.628754 -0.1940.460237-0.1140.633805 -0.1920.463934-0.1120.638895 -0.190.46766-0.110.644027 -0.1880.471416-0.1080.6492 -0.1860.475203-0.1060.654414 -0.1840.47902-0.1040.659671 -0.1820.482867-0.1020.664969 -0.180.486746-0.10.67031 -0.1780.490655-0.0980.675694 -0.1760.494596-0.0960.681121 -0.1740.498569-0.0940.686592 -0.1720.502573-0.0920.692107 -0.170.50661-0.090.697666 -0.1680.510679-0.0880.70327 -0.1660.514781-0.0860.708919 -0.1640.518916-0.0840.714613 -0.1620.523084-0.0820.720352 -0.160.527285-0.080.726138 -0.1580.53152-0.0780.731971 -0.1560.535789-0.0760.73785 -0.1540.540093-0.0740.743776 -0.1520.544431-0.0720.749751-0.070.7557730.01 1.040796 -0.0680.7618430.012 1.049156-0.0660.7679620.014 1.057583 -0.0640.7741310.016 1.066077 -0.0620.7803480.018 1.07464 -0.060.7866160.02 1.083272 -0.0580.7929340.022 1.091973 -0.0560.7993030.024 1.100743 -0.0540.8057230.026 1.109585 -0.0520.8121950.028 1.118497 -0.050.8187190.03 1.127481 -0.0480.8252950.032 1.136537 -0.0460.8319240.034 1.145666 -0.0440.8386060.036 1.154868 -0.0420.8453410.038 1.164144 -0.040.8521310.04 1.173494 -0.0380.8589760.042 1.18292 -0.0360.8658750.044 1.192421 -0.0340.872830.046 1.201999 -0.0320.879840.048 1.211654 -0.030.8869070.05 1.221386 -0.0280.8940310.052 1.231196 -0.0260.9012120.054 1.241085 -0.0240.9084510.056 1.251054 -0.0220.9157480.058 1.261102 -0.020.9231030.06 1.271232 -0.0180.9305170.062 1.281442 -0.0160.9379910.064 1.291735 -0.0140.9455250.066 1.30211 -0.0120.953120.068 1.312569 -0.010.9607750.07 1.323112 -0.0080.9684930.072 1.333739 -0.0060.9762720.074 1.344452 -0.0040.9841130.076 1.355251 -0.0020.9920180.078 1.366136 00.9999860.08 1.377109 0.002 1.0080180.082 1.38817 0.004 1.0161140.084 1.39932 0.006 1.0242760.086 1.41056 0.008 1.0325030.088 1.4218890.09 1.433310.17 1.973853 0.092 1.4448230.172 1.9897080.094 1.4564280.174 2.005689 0.096 1.4681260.176 2.021799 0.098 1.4799180.178 2.038039 0.1 1.4918050.18 2.054408 0.102 1.5037870.182 2.07091 0.104 1.5158660.184 2.087543 0.106 1.5280410.186 2.104311 0.108 1.5403150.188 2.121213 0.11 1.5526870.19 2.138251 0.112 1.5651580.192 2.155425 0.114 1.577730.194 2.172738 0.116 1.5904020.196 2.19019 0.118 1.6031760.198 2.207781 0.12 1.6160530.2 2.225515 0.122 1.6290340.202 2.24339 0.124 1.6421180.204 2.261409 0.126 1.6553080.206 2.279573 0.128 1.6686040.208 2.297883 0.13 1.6820060.21 2.31634 0.132 1.6955160.212 2.334945 0.134 1.7091350.214 2.3537 0.136 1.7228630.216 2.372605 0.138 1.7367010.218 2.391662 0.14 1.750650.22 2.410872 0.142 1.7647120.222 2.430236 0.144 1.7788860.224 2.449756 0.146 1.7931740.226 2.469433 0.148 1.8075770.228 2.489268 0.15 1.8220960.23 2.509262 0.152 1.8367310.232 2.529417 0.154 1.8514840.234 2.549733 0.156 1.8663550.236 2.570213 0.158 1.8813460.238 2.590857 0.16 1.8964570.24 2.611667 0.162 1.911690.242 2.632645 0.164 1.9270450.244 2.65379 0.166 1.9425230.246 2.675106 0.168 1.9581260.248 2.6965930.25 2.7182520.33 3.743385 0.252 2.7400850.332 3.7734530.254 2.7620940.334 3.803761 0.256 2.7842790.336 3.834314 0.258 2.8066430.338 3.865111 0.26 2.8291860.34 3.896156 0.262 2.8519110.342 3.927451 0.264 2.8748180.344 3.958996 0.266 2.8979090.346 3.990796 0.268 2.9211850.348 4.02285 0.27 2.9446480.35 4.055162 0.272 2.96830.352 4.087734 0.274 2.9921420.354 4.120567 0.276 3.0161750.356 4.153664 0.278 3.0404010.358 4.187026 0.28 3.0648220.36 4.220657 0.282 3.0894390.362 4.254558 0.284 3.1142540.364 4.288731 0.286 3.1392680.366 4.323179 0.288 3.1644830.368 4.357903 0.29 3.18990.37 4.392906 0.292 3.2155220.372 4.42819 0.294 3.2413490.374 4.463758 0.296 3.2673840.376 4.499612 0.298 3.2936280.378 4.535753 0.3 3.3200830.38 4.572185 0.302 3.346750.382 4.608909 0.304 3.3736320.384 4.645928 0.306 3.4007290.386 4.683245 0.308 3.4280440.388 4.720861 0.31 3.4555790.39 4.75878 0.312 3.4833350.392 4.797003 0.314 3.5113130.394 4.835533 0.316 3.5395160.396 4.874373 0.318 3.5679460.398 4.913524 0.32 3.5966040.4 4.95299 0.322 3.6254930.402 4.992773 0.324 3.6546130.404 5.032876 0.326 3.6839670.406 5.0733 0.328 3.7135570.408 5.114050.41 5.1551260.497.099276 0.412 5.1965330.4927.1562980.414 5.2382720.4947.213778 0.416 5.2803460.4967.27172 0.418 5.3227590.4987.330127 0.42 5.3655120.57.389004 0.422 5.4086080.5027.448353 0.424 5.4520510.5047.508179 0.426 5.4958420.5067.568486 0.428 5.5399850.5087.629277 0.43 5.5844830.517.690556 0.432 5.6293380.5127.752327 0.434 5.6745540.5147.814595 0.436 5.7201330.5167.877362 0.438 5.7660770.5187.940634 0.44 5.8123910.528.004414 0.442 5.8590770.5228.068707 0.444 5.9061380.5248.133516 0.446 5.9535770.5268.198845 0.448 6.0013960.5288.264699 0.45 6.04960.538.331082 0.452 6.0981910.5328.397998 0.454 6.1471730.5348.465452 0.456 6.1965480.5368.533447 0.458 6.2463190.5388.601989 0.46 6.296490.548.671081 0.462 6.3470640.5428.740728 0.464 6.3980450.5448.810935 0.466 6.4494340.5468.881705 0.468 6.5012370.5488.953044 0.47 6.5534560.559.024956 0.472 6.6060940.5529.097445 0.474 6.6591550.5549.170517 0.476 6.7126420.5569.244175 0.478 6.7665580.5589.318426 0.48 6.8209080.569.393272 0.482 6.8756950.5629.46872 0.484 6.9309210.5649.544774 0.486 6.9865910.5669.621439 0.4887.0427080.5689.6987190.579.776620.6513.46367 0.5729.8551470.65213.571810.5749.9343050.65413.68082 0.57610.01410.65613.79071 0.57810.094530.65813.90147 0.5810.175610.6614.01313 0.58210.257340.66214.12569 0.58410.339730.66414.23915 0.58610.422780.66614.35352 0.58810.50650.66814.46881 0.5910.590890.6714.58502 0.59210.675960.67214.70217 0.59410.761710.67414.82026 0.59610.848150.67614.9393 0.59810.935280.67815.05929 0.611.023110.6815.18025 0.60211.111650.68215.30218 0.60411.20090.68415.42509 0.60611.290870.68615.54898 0.60811.381560.68815.67387 0.6111.472980.6915.79977 0.61211.565130.69215.92667 0.61411.658020.69416.0546 0.61611.751660.69616.18355 0.61811.846050.69816.31354 0.6211.94120.716.44457 0.62212.037110.70216.57665 0.62412.133790.70416.7098 0.62612.231250.70616.84401 0.62812.32950.70816.97931 0.6312.428530.7117.11569 0.63212.528360.71217.25316 0.63412.628990.71417.39174 0.63612.730420.71617.53143 0.63812.832680.71817.67225 0.6412.935750.7217.81419 0.64213.039650.72217.95728 0.64413.144390.72418.10151 0.64613.249960.72618.24691 0.64813.356390.72818.393470.7318.541210.8125.53363 0.73218.690130.81225.738720.73418.840250.81425.94545 0.73618.991580.81626.15385 0.73819.144120.81826.36392 0.7419.297890.8226.57568 0.74219.452890.82226.78914 0.74419.609140.82427.00431 0.74619.766640.82627.22121 0.74819.925410.82827.43985 0.7520.085450.8327.66025 0.75220.246780.83227.88242 0.75420.409410.83428.10638 0.75620.573340.83628.33213 0.75820.738580.83828.5597 0.7620.905160.8428.78909 0.76221.073070.84229.02033 0.76421.242330.84429.25342 0.76621.412950.84629.48839 0.76821.584940.84829.72524 0.7721.758310.8529.964 0.77221.933080.85230.20467 0.77422.109250.85430.44728 0.77622.286830.85630.69184 0.77822.465840.85830.93836 0.7822.646290.8631.18686 0.78222.828190.86231.43735 0.78423.011550.86431.68986 0.78623.196380.86631.9444 0.78823.382690.86832.20098 0.7923.570510.8732.45962 0.79223.759830.87232.72034 0.79423.950670.87432.98315 0.79624.143040.87633.24807 0.79824.336960.87833.51513 0.824.532440.8833.78432 0.80224.729490.88234.05568 0.80424.928110.88434.32922 0.80625.128340.88634.60496 0.80825.330170.88834.882910.8935.163090.94643.99154 0.89235.445520.94844.344880.89435.730220.9544.701070.89636.017210.95245.060110.89836.306510.95445.422040.936.598120.95645.786870.90236.892080.95846.154630.90437.188410.9646.525350.90637.487110.96246.899050.90837.788210.96447.275750.9138.091730.96647.655470.91238.397680.96848.038240.91438.70610.9748.424090.91639.016990.97248.813040.91839.330380.97449.205110.9239.646280.97649.600330.92239.964720.97849.998720.92440.285720.9850.400320.92640.60930.98250.805140.92840.935480.98451.213210.9341.264280.98651.624560.93241.595720.98852.039210.93441.929820.9952.45720.93642.26660.99252.878540.93842.606090.99453.303270.9442.948310.99653.73140.94243.293270.99854.162980.94443.64101154.59802复化Simpson数据：-1 0.018319929 -0.34 0.256658088 0.32 3.596641805 -0.98 0.0198445 -0.32 0.278035042 0.34 3.896195298-0.96 0.021494322 -0.3 0.301192133 0.36 4.220697765-0.94 0.023283225 -0.28 0.326278124 0.38 4.572227037-0.92 0.025220379 -0.26 0.353453177 0.4 4.95303418-0.9 0.027320224 -0.24 0.382891765 0.42 5.365557596-0.88 0.029594431 -0.22 0.41478194 0.44 5.812438891-0.86 0.032059069 -0.16 0.527292277 0.54 8.671138204-0.84 0.034728638 -0.14 0.571209036 0.56 9.39333156-0.82 0.037621263 -0.12 0.61878367 0.58 10.17567433-0.8 0.040754615 -0.1 0.670320427 0.6 11.02317608-0.78 0.044149394 -0.08 0.726149698 0.62 11.94126383-0.76 0.047826844 -0.06 0.78662861 0.64 12.93581634-0.74 0.051810827 -0.04 0.85214479 0.66 14.01320231-0.72 0.056126648 -0.02 0.92311742 0.68 15.1803205-0.7 0.060802006 0 1.0000013 0.7 16.44464467 -0.68 0.065866854 0.02 1.083288424 0.72 17.81427057 -0.66 0.071353499 0.04 1.173512427 0.74 19.29796874 -0.64 0.077297255 0.06 1.271250748 0.76 20.90523965 -0.62 0.083735917 0.08 1.377129533 0.78 22.64637562 -0.6 0.090711017 0.1 1.491826493 0.8 24.53252554 -0.58 0.098266855 0.12 1.616076341 0.82 26.57576756 -0.56 0.106452202 0.14 1.750674449 0.84 28.78918506 -0.54 0.11531904 0.16 1.896482943 0.86 31.18695183 -0.52 0.12492459 0.18 2.054435268 0.88 33.78442141 -0.5 0.135329888 0.2 2.225543071 0.9 36.59822683 -0.48 0.14660204 0.22 2.410901825 0.92 39.64638571 -0.46 0.158812728 0.24 2.611698647 0.94 42.94841704 -0.44 0.17204064 0.26 2.829219145 0.96 46.52546475 -0.42 0.18636997 0.28 3.064856356 0.98 50.40043451 -0.4 0.201892977 0.3 3.320119013 1 54.59813904 -0.38 0.218708553 0.46 6.296539601-0.36 0.236924875 0.48 6.820959636-0.2 0.449328351 0.5 7.389057081-0.18 0.486751777 0.52 8.0044696750102030405060四、讨论①在满足相同精度要求的情况下复化梯形积分法比复化Simpson 积分法计算所需节点数多，计算量大。

北京航空航天大学2020届研究生《数值分析》实验作业第九题院系：xx学院学号：姓名：2020年11月Q9:方程组A.4一、 算法设计方案（一）总体思路1.题目要求∑∑===k i kj s r rsy x cy x p 00),(对f(x, y) 进行拟合，可选用乘积型最小二乘拟合。

),(i i y x 与),(i i y x f 的数表由方程组与表A-1得到。

2.),(**j i y x f 与1使用相同方法求得，),(**j i y x p 由计算得出的p(x,y)直接带入),(**j i y x 求得。

1. ),(i i y x 与),(i i y x f 的数表的获得对区域D ={ （x,y)|1≤x ≤1.24,1.0≤y ≤1.16}上的f (x , y )值可通过xi=1+0.008i ，yj=1+0.008j ，得到),(i i y x 共31×21组。

将每组带入A4方程组，即可获得五个二元函数组，通过简单牛顿迭代法求解这五个二元数组可获得z1~z5有关x,y 的表达式。

再将),(i i y x 分别带入z1~z5表达式即可获得f(x,y)值。

2.乘积型最小二乘曲面拟合2.1使用乘积型最小二乘拟合，根据k 值不用，有基函数矩阵如下：⎪⎪⎪⎭⎫ ⎝⎛=k i i k x x x x B 0000 ， ⎪⎪⎪⎭⎫ ⎝⎛=k j jk y y y y G 0000数表矩阵如下：⎪⎪⎪⎭⎫⎝⎛=),(),(),(),(0000j i i j y x f y x f y x f y x f U记C=[rs c ]，则系数rs c 的表达式矩阵为：11-)(-=G G UG B B B C T TT ）（通过求解如下线性方程，即可得到系数矩阵C 。

UG B G G C B B T T T =)()(2.2计算),(),,(****j i j i y x p y x f (i =1,2,…,31 ; j =1,2,…,21) 的值),(**j i y x f 的计算与),(j i y x f 相同。

《数值分析》作业三院系：机械学院学号：SY1307145姓名：龙安林2013年11 月24 日1. 算法设计1) 开始；2) 计算数组[][]0.08,0.050.5,0,1,2,,10;0,1,2,,20x i i y j j i j ==+=⋯=⋯（）； 3) 将点[][],0,1,2,,10;0,1,2,,20x i y j i j =⋯=⋯（）,（）带入非线性方程组： 0.5cos 2.670.5sin 1.070.5cos 3.740.5sin 0.79t u v w x t u v w y t u v w x t u v w y +++-=⎧⎪+++-=⎪⎨+++-=⎪⎪+++-=⎩ 得出相应的点,t u （）； 4) 选择拉格朗日插值法，将,t u （）作为中间变量，在题目所给出的二维数表中进行二次代数插值，得到[][],)(z f x i y j =；5) 输出数表：[][][][]()()0,1,2,,10;0,1,2,,20,,,x i y j f x i y j i j =⋯=⋯； 6) 令k=0；7) 以()()(),,,0,1,r r r s x x y y r s ϕψ===…,k 为拟合基函数，将上述数表作为拟合条件，对于给定的k 值，得到矩阵B 、G 、U ；8) 令-1-1(),()T T T A B B B U C AG G G ==，用选主元的LU 分解法分别计算矩阵A 和C 的各列，最后得到系数矩阵C ；9) 以公式：()()()00,k ki j rs r i s j s r p x y C x y ϕψ===∑∑计算每个点的拟合值；10) 利用公式：()()()2102000,,i j i j i j f x y p x y σ===-∑∑计算拟合误差，当σ≤10-7时，循环结束，否则k=k+1，转（6）；11) 令[][]()**0.10.50.2 1,2,81,2,5x i i y j j i j ==+=⋯=⋯；，；，；12) 计算()()()******,,,,,i j i j i jf x y p x y delta x y ，输出数表，观察逼近效果； 13) 结束。

数值分析大作业数值分析大作业姓名：黄晨晨学号：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：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=时，计算的结果如下：此结果即为上文中的正确计算结果。

数值分析—计算实习作业一学院：机械工程学院专业：材料加工工程姓名：暴一品学号：SY12071342012-10-29一、算法设计方案观察矩阵A ，结构为带状，且与主对角线相邻的两个带的值b 和c 都是常数。

从而可以用带原点平移的幂法或反幂法计算λ1和λ501。

所以算法的设计方案如下：1.求按模最大的特征值，并记为max_eigenvalue ，算法如下所示⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=======------≤≤-),2,1()sgn(),,(/max ),,()(1)()(11)1(11)1(1)1()0()0(10ΛΛΛk h h h h Ay u h u y h h h h u k r k r k Tk nk k kk r k k k j nj k rTn β任取非零向量2.平移矩阵得到A ’=A-max_eigenvalueI ，再次用幂法，这次求出的A ’的按模大的特征值pymax_eigenvalue 就是与步骤1求出的特征值相差最大的特征值。

即两者一个为最大的特征值，另一个为最小的特征值。

3.根据max_eigenvalue 和pymax_eigenvalue 的正负性，直接确定λ1，和λ501。

4.对原矩阵A 用反幂法，求出其按模最小的特征值，记为s_eigenvalue ，此即λs 。

⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=====∈--------),2,1(/111111110Λk u y y Au u y u u R u k T k k k k k k k k Tk k n βηη任取非零向量在反幂法的求解过程中，每迭代一次都要求满足解线性方程组Auk=yk-1。

本题中矩阵A 上半带宽为2，下半带宽也为2 。

故选择采用三角分解法求解方程组：先将原矩阵改写成5行501列的矩阵C （不存储A 的0元素） A 的带内元素aij=c 中的元素ci-j+3。

再对C 矩阵做LU 分解。

对于k=1,2，…，n ，执行∑---=+-+-+-+--=1)2,2,1max(,3,3,3,3:k j k t jj t t t k j j k j j k ccc c [j=k,k+1,…,min(k+2,n)]kk s k r i t k k t t t i k k i k k i c ccc c ,31),,1max(,3,3,3,3/)(:∑---=+-+-+-+--=[i=k+1,k+2,…,min(k+2,n);k<n]求解Lx=b ，Uuk=x （数组b 先是存放原方程组右端向量yk-1，后来存放中间向量x ）∑--=+--=1),1max(,3:i r i t tt t i i i bcb b （i=2,3,…,n ）nn kn c b u ,3/:=in i i t kt tt i i ki c u cb u ,3),2min(1,3/)(:∑++=+--= (i=n-1,n-2, (1)5.对k=1，2，……39执行：先根据题中给出的公式算出μk ，再将矩阵平移A ”=A-μk ，对矩阵A ”运用反幂法（线性方程组的解法同上），就可以求出与μk 最接近的特征值λik ，保存在数组py_eigenvalue 中。

Hilbert 矩阵病态问题研究 （数值分析第一次大作业） 姓名：** 学号：** 班级：**1）Hilbert 矩阵的阶数n 与ln(())n cond H 的关系猜想：ln(())n cond H 与n 呈线性关系，其中()n cond H 按2范数计算。

绘制ln(())n cond H n 曲线。

分别取11050500n ≤≤、、，得到ln(())n cond H n 曲线如图1-1、图1-2及图1-3所示。

程序详见附录1。

图1-1. 110n ≤≤由图1-1可知，110n ≤≤，ln(())n cond H 是n 的线性函数，猜想正确。

图1-2. 150n ≤≤由图1-2知，当15n >时，ln(())n cond H 与n 之间的线性关系已经不存在，而且ln(())n cond H 的值大致在(40,50)内间波动，猜想与实际不完全相符。

图1-3. 1500n ≤≤图1-3进一步说明了ln(())n cond H 与n 之间的变化关系：当n 小于某一值（设该值为k ）时，ln(())n cond H 是n 的线性函数，而当n 大于k 时，随着n 的增大，ln(())n cond H 与n 间的线性关系不再成立，且其值在某一区间内波动。

为进一步确定k 的大小，绘制114n ≤≤时的曲线，如图1-4所示，可知k 的取值应为13。

图1-4. 114n ≤≤2）由n H 至ˆnH 的预处理 绘制ˆln(()/())n n cond H cond H n 曲线。

其中11ˆn nH D H D --=，D 为由n H 的对角元素开方构成的对角矩阵。

条件数按2范数计算。

程序详见附录2。

分别取11350500n ≤≤、、，得到如图2-1、图2-2和图2-3所示曲线。

由曲线图像可知：当Hilbert 矩阵的阶数12n ≤时，ˆln(()/())n ncond H cond H 随n 增大而逐渐减小，而n 继续增大时，ˆln(()/())n n cond H cond H 的取值将在区间(-7,4)内波动，且主要集中在(0,-3)区间内。

数值分析大作业一、算法设计方案1、矩阵初始化矩阵[]501501⨯=ij a A 的下半带宽r=2,上半带宽s=2，设置矩阵[][]5011++s r C ，在矩阵C 中检索矩阵A 中的带内元素ij a 的方法是：j s j i ij c a ,1++-=。

这样所需要的存储单元数大大减少，从而极大提高了运算效率。

2、利用幂法求出5011λλ，幂法迭代格式：0111111nk k k k kk T k k k u R y u u Ay y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止。

首先对于矩阵A 利用幂法迭代求出一个λ，然后求出矩阵B ，其中I A B λ-=(I 为单位矩阵），对矩阵B 进行幂法迭代，求出λ'，之后令λλλ+'=''，比较的大小与λλ''，大者为501λ，小者为1λ。

3、利用反幂法求出ik s λλ,反幂法迭代格式：0111111nk k k k kk T k k k u R y u Au y y u ηηβ------⎧∈⎪⎪=⎪=⎨⎪=⎪⎪=⎩非零向量 当1210/-≤-k k βββ时，迭代终止，1s k λβ=。

每迭代一次都要求解一次线性方程组1-=k k y Au ，求解过程为：（1）作分解LU A =对于n k ,...,2,1=执行[][]s k n r k k k i c c c c c n s k k k j c cc c k s ks k t k s k r i t t s t i k s k i k s k i js j t k s j r k t t s t k j s j k j s j k <+++=-=++=-=+++----=++-++-++-++----=++-++-++-∑∑);,min(,...,2,1/)(:),min(,...,1,:,1,11),,1max(,1,1,1,11),,1max(,1,1,1（2）求解y Ux b Ly ==,（数组b 先是存放原方程组右端向量，后来存放中间向量y))1,...,2,1(/)(:/:),...,3,2(:,1),min(1.1.11),1max(,1--=-===-=+++-++-+--=++-∑∑n n i c x c b x c b x n i b c b b i s t n s i i t t s t i i i ns n n ti r i t t s t i i i使用反幂法，直接可以求得矩阵按模最小的特征值s λ。

数值分析上机作业第 1 章1.1计算积分，n=9。

（要求计算结果具有6位有效数字）程序：n=1:19;I=zeros(1,19);I(19)=1/2*((exp(-1)/20)+(1/20));I(18)=1/2*((exp(-1)/19)+(1/19));for i=2:10I(19-i)=1/(20-i)*(1-I(20-i));endformat longdisp(I(1:19))结果截图及分析：在MATLAB中运行以上代码，得到结果如下图所示：当计算到数列的第10项时，所得的结果即为n=9时的准确积分值。

取6位有效数字可得.1.2分别将区间[-10.10]分为100,200,400等份，利用mesh或surf命令画出二元函数z=的三维图形。

程序：>> x = -10:0.1:10;y = -10:0.1:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.1')>> x = -10:0.2:10;y = -10:0.2:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长 0.2')>>x = -10:0.05:10;y = -10:0.05:10;[X,Y] = meshgrid(x,y);Z = exp(-abs(X))+cos(X+Y)+1./(X.^2+Y.^2+1);subplot(2,2,1);mesh(X,Y,Z);title('步长0.05')结果截图及分析：由图可知，步长越小时，绘得的图形越精确。

实验报告课程名称：数值分析实验项目：曲线拟合/数值积分专业班级：姓名：学号：实验室号：实验组号：实验时间：20.10.24 批阅时间：指导教师：成绩：工业大学实验报告（适用计算机程序设计类）专业班级：学号：姓名：实验名称：曲线拟合与函数插值附件A 工业大学实验报告（适用计算机程序设计类）专业班级：学号：姓名：实验步骤或程序：附录一：1.利用二次，三次，四次多项式进行拟合：1.1 MATLAB代码如下：clear;clc;close allt=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24];y=[14 13 13 13 13 14 15 17 19 21 22 24 27 30 31 30 28 26 24 23 21 19 17 16 15];%输入数据hold on[p2 s2]=polyfit(t,y,2);%对于上面的数据进行2次多项式拟合,其中s2包括R（系数矩阵的QR分解的上三角阵），%df（自由度），normr（拟合误差平方和的算术平方根）。

y2=polyval(p2,t);%返回多项式拟合曲线在t处的值[p3 s3]=polyfit(t,y,3);y3=polyval(p3,t);[p4 s4]=polyfit(t,y,4);y4=polyval(p4,t);plot(t,y,'ro')%画图plot(t,y2,'g-')plot(t,y3,'m^-')plot(t,y4,'bs-')xlabel('t')ylabel('y')legend('原始数据','2次多项式拟合','3次多项式拟合','4次多项式拟合')1.2 二次，三次，四次多项式拟合的结果分别如下：（1）总的拟合结果在工作区的显示如下：（2）其次二次多项式拟合的结果为：（3）其中三次多项式拟合的结果：（4）其中四次多项式拟合的结果为：1.3 拟合的图像为：1.4 拟合的多项式为：根据工作区得出的数据列出最后的拟合多项式为：(1)y=7.416+2.594t-0.094t^2(2)y=12.251-0.102t+0.193t^2-0.008t^3(3)y=15.604-3.526t+0.866t^2-0.052t^3+0.0009t^42.形如2()()b t c y t ae--=的函数，其中,,a b c 为待定常数。

数值分析大作业数值分析大作业学号：*********专业：机械工程学生姓名：***2014年10月摘要：在自然科学与工程技术中，很多问题的解决常常归结为求解线性方程组Ax=b 。

随着计算机的发展，利用计算机这个强有力的计算工具去求解线性方程组是一个非常实用的问题。

在求解大型线性方程组时，直接法在多次消元，回代的过程中，四则运算的误差累计与传播无法控制，致使计算结果的精度就无法保证，特别是求解大型稀松矩阵时，还要对系数矩阵进行分解。

而迭代法相对于直接法而言，具有保持迭代矩阵不变的特点，计算程序一般也比较简单，且对于许多问题收敛速度比较快。

比较常用的迭代法有雅克比迭代法、高斯一塞德尔迭代法和逐次超松弛迭代法等，本次研究目的是通过求解一个线性方程组来比较它们的迭代效果，验证一些已有的结论。

1．数学原理1.1雅可比迭代法将线性方程组的系数矩n n ij R a A ?∈=)(分解为A=D+L+U ，其中D 是由A 的主对角元素构成的对角矩阵，L 是由A 的严格下三角部分构成的严格下三角矩阵， U 是由A 的严格上三角部分构成的严格上三角矩阵，即,2211=nn a a a D.0000,0000,1223113121,21323121==--n n n n n n n n a a a a a a U a a a a a a L若系数矩阵A 的对角元素),,2,1(0n i a ii =≠，则矩阵D 非奇异，取M=D ，N=-（L+U)，则J J g x G b D x U L D x +=++-=--11)(，因而，构造的迭代法为：.),(,11)()1(b D g U L D G g x G x J J J k J k --+=+-=+=1.2高斯-赛得尔迭代法将线性方程组的系数矩n n ij R a A ?∈=)(分解为A=D+L+U 。

若系数矩阵A 的对角元素不等于0，则矩阵D 非奇异，取M=L+D ，N=-U ，则()()G G g x G b D L Ux D L x +=+++-=--11因而，构造的迭代法为：()().,,11)()1(b D L g U D L G g x G x G G G k G k --++=+-=+=1.3逐次超松弛迭代法线性方程组的系数矩n n ij R a A ?∈=)(分解为A=D+L+U 。