数值分析大作业

数值分析大作业学号:*********专业:机械工程学生姓名:***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 。

取()()[],11,1U D N L D M ωωωωω--=+=则()()[]()S S g x G b L D x U D L D x +=++--+=--111ωωωωω因而,构造的迭代法为:()()[]().,1,11)()1(b L D g U D L D G g x G x S G S k S k --++=--+=+=ωωωωω2.程序设计2.1雅可比迭代法求解雅可比迭代法MATLAB 程序如下: %majacobi.mfunction x=majacobi(A,b,x0,ep,N)%用途:用Jacobi 迭代法解线性方程组Ax=b%格式;x=majacobi(A,b,x0,ep,N) A 为系数矩阵,b 为右端向量, %x0为初始向量(默认零向量),ep 为精度(默认1e-6), %N 为最大迭代次数(默认500次),x 返回近似解向量 n=length(b);if nargin<5,N=500;end if nargin<4,ep=1e-6;endif nargin<3,x0=zeros(n,1);end x=zeros(n,1);k=0; while k<N for i=1:nx(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); endif norm(x-x0,inf)<ep,break; end x0=x ;k=k+1; endif k==N,Warning('100');end disp(['k=',num2str(k)])在MATLAB 命令窗口执行创建的m 文件 >> edit majacobi>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A =10.9000 1.2000 2.1000 0.90001.2000 11.2000 1.50002.50002.1000 1.5000 9.8000 1.30000.9000 2.5000 1.3000 12.3000>> b=[-7.0,5.3,10.3,24.6]';>> x=majacobi(A,b)k=17x =-0.99860.00711.00321.96562.2高斯-赛得尔迭代法求解高斯-赛得尔迭代法MATLAB程序如下:%maseidel.mfunction x=maseidel (A,b,x0,ep,N)%用途:用Gauss-Seidel迭代法解线性方程组Ax=b%格式:x=maseidel (A,b,x0,ep,N) A为系数矩阵,b为右端向量,%x0为初始向量(默认零向量),ep为精度(默认1e-6),%N为最大迭代次数(默认500次),x返回近似解向量n=length(b);if nargin<5,N=500;endif nargin<4,ep=1e-6;endif nargin<3,x0=zeros(n,1);endx=zeros(n,1);k=0;while k<Nfor i=1:nif i==1x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);else if i==nx(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);elsex(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endendendif norm(x-x0,inf)<ep,break;endx0=x;k=k+1;endif k==N,Warning('100');enddisp(['k=',num2str(k)])在MATLAB命令窗口执行创建的m文件:>> edit maseidel>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A =10.9000 1.2000 2.1000 0.90001.2000 11.2000 1.50002.50002.1000 1.5000 9.8000 1.30000.9000 2.5000 1.3000 12.3000>> b=[-7.0,5.3,10.3,24.6]';>> x=maseidel(A,b)k=7x =-0.99860.00711.00321.96562.3逐次超松弛迭代法求解逐次超松弛迭代法MATLAB程序如下:%masor.mfunction x=masor(A,b,omega,x0,ep,N)%用途:用SOR迭代法解线性方程组Ax=b%格式:x=maseidel (A,b,x0,ep,N) A为系数矩阵,b为右端向量,%omega为松弛因子(默认1.5),x0为初始向量(默认零向量),ep为精度(默认1e-6),%N为最大迭代次数(默认500次),x返回近似解向量n=length(b);if nargin<6,N=500;endif nargin<5,ep=1e-6;endif nargin<4,x0=zeros(n,1);endif nargin<3,omega=1.5;endx=zeros(n,1);k=0;while k<Nfor i=1:nif i==1x1(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);else if i==nx1(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);elsex1(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);endendx(i)=(1-omega)*x0(i)+omega*x1(i);endif norm(x0-x,inf)<ep,break;endk=k+1;x0=x;endif k==N,Warning('100');enddisp(['k=',num2str(k)])取不同的松弛因子进行试验,在MATLAB命令窗口执行创建的m文件:>> edit masor>> A=[10.9,1.2,2.1,0.9;1.2,11.2,1.5,2.5;2.1,1.5,9.8,1.3;0.9,2.5,1.3,12.3]A =10.9000 1.2000 2.1000 0.90001.2000 11.2000 1.50002.50002.1000 1.5000 9.8000 1.30000.9000 2.5000 1.3000 12.3000>> b=[-7.0,5.3,10.3,24.6]';>> x=masor(A,b,0.4)k=33x =-0.99860.00711.00321.9656取其它的松弛因子得:取松弛因子为0.8 >> x=masor(A,b,0.8) 得到k=12;x = -0.9986,0.0071,1.0032,1.9656。

取松弛因子为1.2 >> x=masor(A,b,1.2) 得到k=11;x =-0.9986,0.0071,1.0032,1.9656。

取松弛因子为1.6 >> x=masor(A,b,1.6) 得到k=33;x =-0.9986,0.0071,1.0032,1.9656。

三.结果分析和讨论由程序解得的结果可知:雅可比迭代法在第17次迭代达到收敛,高斯-塞德尔迭代法在第7次迭代达到收敛。

从此例可以看出,高斯-塞德尔迭代法比雅可比迭代法收敛速度快(达到同样的精度所需迭代次数少)。

由雅可比迭代公式可知,在迭代的每一步计算过程中使用的是旧的分量来计算,没有用到迭代出的新的分量。

合集下载

北航数值分析大作业 第一题 幂法与反幂法

北航数值分析大作业 第一题 幂法与反幂法

数 值 分 析(B ) 大 作 业(一)姓名: 学号: 电话:1、算法设计:①求1λ、501λ和s λ的值:s λ:s λ表示矩阵的按模最小特征值,为求得s λ直接对待求矩阵A 应用反幂法即可。

1λ、501λ:若矩阵A 的特征值满足关系 1n λλ<<且1n λλ≠,要求1λ、及501λ时,可按如下方法求解: a . 对矩阵A 用幂法,求得按模最大的特征值1m λ。

b . 按平移量1m λ对矩阵A 进行原点平移得矩阵1m BA I λ=+,对矩阵B 用反幂法求得B 的按模最小特征值2m λ。

c . 321m m m λλλ=-则:113min(,)m m λλλ=,13max(,)n m m λλλ=即为所求。

②求和A 的与数5011140k k λλμλ-=+最接近的特征值ik λ(k=0,1,…39):求矩阵A 的特征值中与P 最接近的特征值的大小,采用原点平移的方法:先求矩阵 B=A-PI 对应的按模最小特征值k β,则k β+P 即为矩阵A 与P 最接近的特征值。

在本次计算实习中则是先求平移矩阵k B A I μ=-,对该矩阵应用反幂法求得s λ,则与k μ最接近的A 的特征值为:s P λ+重复以上过程39次即可求得ik λ(k=0,1,…39)的值。

③求A 的(谱范数)条件数2cond()A 和行列式det A :在(1)中用反幂法求矩阵A 的按模最小特征值时,要用到Doolittle 分解方法,在Doolittle 分解完成后得到的两个矩阵分别为L 和U ,则A 的行列式可由U 阵求出,即:det(A)=det(U)。

求得det(A)不为0,因此A 为非奇异的实对称矩阵,则: max 2()scond A λλ=,max λ和s λ分别为模最大特征值与模最小特征值。

2、程序源代码:#include "Stdio.h"#include "Conio.h"#include "math.h"//****************************************************************************// // 在存储带状矩阵时,下面的几个量在程序中反复用到,为方便编程故把它们定义成宏.// // M :转换后的矩阵的行数,M=R+S+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];}}}}三、程序结果。

数值分析大作业

数值分析大作业

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

北航数值分析全部三次大作业

北航数值分析全部三次大作业

北航数值分析全部三次大作业第一次大作业是关于解线性方程组的数值方法。

我们被要求实现各种常用的线性方程组求解算法,例如高斯消元法、LU分解法和迭代法等。

我首先学习了这些算法的原理和实现方法,并借助Python编程语言编写了这些算法的代码。

在实验中,我们使用了不同规模和条件的线性方程组进行测试,并比较了不同算法的性能和精度。

通过这个作业,我深入了解了线性方程组求解的原理和方法,提高了我的编程和数值计算能力。

第二次大作业是关于数值积分的方法。

数值积分是数值分析中的重要内容,它可以用于计算曲线的长度、函数的面积以及求解微分方程等问题。

在这个作业中,我们需要实现不同的数值积分算法,例如矩形法、梯形法和辛普森法等。

我学习了这些算法的原理和实现方法,并使用Python编写了它们的代码。

在实验中,我们计算了不同函数的积分值,并对比了不同算法的精度和效率。

通过这个作业,我深入了解了数值积分的原理和方法,提高了我的编程和数学建模能力。

第三次大作业是关于常微分方程的数值解法。

常微分方程是数值分析中的核心内容之一,它可以用于描述众多物理、化学和生物现象。

在这个作业中,我们需要实现不同的常微分方程求解算法,例如欧拉法、龙格-库塔法和Adams法等。

我学习了这些算法的原理和实现方法,并使用Python编写了它们的代码。

在实验中,我们解决了一些具体的常微分方程问题,并比较了不同算法的精度和效率。

通过这个作业,我深入了解了常微分方程的原理和方法,提高了我的编程和问题求解能力。

总的来说,北航数值分析课程的三次大作业非常有挑战性,但也非常有意义。

通过这些作业,我在数值计算和编程方面得到了很大的提升,也更加深入地了解了数值分析的理论和方法。

虽然这些作业需要大量的时间和精力,但我相信这些努力将会对我未来的学习和工作产生积极的影响。

北航数值分析报告大作业第三题(fortran)

北航数值分析报告大作业第三题(fortran)

北航数值分析报告大作业第三题(fortran)“数值分析“计算实习大作业第三题——SY1415215孔维鹏一、计算说明1、将x i=0.08i,y j=0.5+0.05j分别代入方程组(A.3)得到关于t,u,v,w的的方程组,调用离散牛顿迭代子函数求出与x i,y j对应的t i,u j。

2、调用分片二次代数插值子函数在点(t i,u j)处插值得到z(x i,y j)=f(x i,y j),得到数表(x i,y j,f(x i,y j))。

3、对于k=1,2,3,4?,分别调用最小二乘拟合子函数计算系数矩阵c rs 及误差σ,直到满足精度,即求得最小的k值及系数矩阵c rs。

4、将x i?=0.1i,y j?=0.5+0.2j分别代入方程组(A.3)得到关于t?,u?,v?,w?的的方程组,调用离散牛顿迭代子函数求出与x i?,y j?对应的t i?,u j?,调用分片二次代数插值子函数在点(t i?,u j?)处插值得到z?(x i?,y j?)=f(x i?,y j?);调用步骤3中求得的系数矩阵c rs求得p(x i?,y j?),打印数表(x i?,y j?,f(x i?,y j?),p(x i?,y j?))。

二、源程序(FORTRAN)PROGRAM SY1415215DIMENSIONX(11),Y(21),T(6),U(6),Z(6,6),UX(11,21),TY(11,21),FXY(11,21), C(6,6) DIMENSIONX1(8),Y1(5),FXY1(8,5),PXY1(8,5),UX1(8,5),TY1(8,5)REAL(8) X,Y,T,U,Z,FXY,UX,TY,C,E,X1,Y1,FXY1,PXY1,UX1,TY1OPEN (1,FILE='第三题计算结果.TXT')DO I=1,11X(I)=0.08*(I-1)ENDDODO I=1,21Y(I)=0.5+0.05*(I-1)ENDDO!*****求解非线性方程组,得到z=f(t,u)的函数*******DO I=1,11DO J=1,21CALL DISNEWTON_NONLINEAR(X(I),Y(J),UX(I,J),TY(I,J)) ENDDO ENDDO!*************分片二次插值得到z=f(x,y)***********DO I=1,11DO J=1,21CALL INTERPOLATION(UX(I,J),TY(I,J),FXY(I,J))ENDDO ENDDOWRITE (1,"('数表(x,y,f(x,y)):')")WRITE (1,"(3X,'X',7X,'Y',10X,'F(X,Y)')")DO I=1,11DO J=1,21WRITE(1,'(1X,F5.2,2X,F5.3,2X,E20.13)') X(I),Y(J),FXY(I,J) ENDDOWRITE (1,"('')")ENDDO!***********最小二乘拟合得到P(x,y)**************N=11M=21WRITE (1,'(" ","K和σ分别为:")')DO K=1,20CALL LSFITTING(X,Y,FXY,C,N,M,K,K,E) WRITE (1,'(I3,2X,E20.13)') K-1,EIF(ETA).OR.(A(L,K)==TA)) THENTA=A(L,K)TL=LDO J=K,NT(K,J)=A(K,J)A(K,J)=A(TL,J)A(TL,J)=T(K,J)ENDDOTB(K)=B(K)B(K)=B(TL)B(TL)=TB(K)ENDIF ENDDODO I=K+1,NM(I,K)=A(I,K)/A(K,K)A(I,K)=0DO J=K+1,NA(I,J)=A(I,J)-M(I,K)*A(K,J) ENDDOB(I)=B(I)-M(I,K)*B(K)ENDDOENDDO!回代过程X(N)=B(N)/A(N,N)DO K=N-1,1,-1S=0.0DO J=K+1,NS=S+A(K,J)*X(J)ENDDOX(K)=(B(K)-S)/A(K,K)ENDDORETURNEND!***********求向量的无穷数************ SUBROUTINE NORM(X,N,A) DIMENSION X(N)REAL(8) X,AA=ABS(X(1))DO I=2,NIF(ABS(X(I))>ABS(X(I-1))) THENA=ABS(X(I)) ENDIFENDDORETURNEND!**************分片二次代数插值************** SUBROUTINE INTERPOLATION(U,V,W) PARAMETER (N=6,M=6)DIMENSION X(N),Y(M),Z(M,N),LK(3),LR(3)REAL(8) X,Y,Z,H,TREAL(8) U,V,W,LK,LR !U,V分别为插值点处的坐标,W为插值结果INTEGER R!**********************数据赋值********************** DATA Y/0.0,0.2,0.4,0.6,0.8,1.0/DATA X/0.0,0.4,0.8,1.2,1.6,2.0/DATA Z/-0.5,-0.42,-0.18,0.22,0.78,1.5,&&-0.34,-0.5,-0.5,-0.34,-0.02,0.46,&&0.14,-0.26,-0.5,-0.58,-0.5,-0.26,&&0.94,0.3,-0.18,-0.5,-0.66,-0.66,&&2.06,1.18,0.46,-0.1,-0.5,-0.74,&&3.5,2.38,1.42,0.62,-0.02,-0.5/H=0.4T=0.2!******************计算K,R************************* IF(UX(N-1)-H/2) THENK=N-1ELSEDO I=3,N-2IF((U>X(I)-H/2).AND.(UY(M-1)-T/2) THENR=M-1 ELSEDO J=3,M-2IF((V>Y(J)-T/2).AND.(VN) P=N IF(P>20) P=20IF(Q>M) Q=MIF(Q>20) Q=20XX=0YY=0D1=NAPX(1)=0.0DO I=1,NAPX(1)=APX(1)+X(I)ENDDOAPX(1)=APX(1)/D1DO J=1,MV(1,J)=0.0DO I=1,NV(1,J)=V(1,J)+Z(I,J)ENDDOV(1,J)=V(1,J)/D1ENDDOIF(P>1) THEND2=0.0APX(2)=0.0DO I=1,NG=X(I)-APX(1)D2=D2+G*GAPX(2)=APX(2)+(X(I)-XX)*G*G ENDDO APX(2)=APX(2)/D2BX(2)=D2/D1DO J=1,MV(2,J)=0.0DO I=1,NG=X(I)-APX(1)V(2,J)=V(2,J)+Z(I,J)*G ENDDOV(2,J)=V(2,J)/D2ENDDOD1=D2ENDIFDO K=3,PD2=0.0APX(K)=0.0DO J=1,MV(K,J)=0.0ENDDODO I=1,NG1=1.0G2=X(I)-APX(1)DO J=3,KG=(X(I)-APX(J-1))*G2-BX(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPX(K)=APX(K)+X(I)*G*GDO J=1,M V(K,J)=V(K,J)+Z(I,J)*G ENDDOENDDODO J=1,MV(K,J)=V(K,J)/D2ENDDOAPX(K)=APX(K)/D2BX(K)=D2/D1D1=D2ENDDOD1=MAPY(1)=0.0DO I=1,MAPY(1)=APY(1)+Y(I)ENDDOAPY(1)=APY(1)/D1DO J=1,PU(J,1)=0.0DO I=1,MU(J,1)=U(J,1)+V(J,I) ENDDO U(J,1)=U(J,1)/D1ENDDOIF(Q>1)THEND2=0.0APY(2)=0.0DO I=1,MG=Y(I)-APY(1)D2=D2+G*G APY(2)=APY(2)+(Y(I))*G*G ENDDO APY(2)=APY(2)/D2BY(2)=D2/D1DO J=1,PU(J,2)=0.0DO I=1,MG=Y(I)-APY(1)U(J,2)=U(J,2)+V(J,I)*GENDDOU(J,2)=U(J,2)/D2ENDDOD1=D2ENDIFDO K=3,QD2=0.0APY(K)=0.0DO J=1,PU(J,K)=0.0ENDDODO I=1,MG1=1.0G2=Y(I)-APY(1)DO J=3,KG=(Y(I)-APY(J-1))*G2-BY(J-1)*G1 G1=G2 G2=GENDDOD2=D2+G*GAPY(K)=APY(K)+Y(I)*G*G DO J=1,PU(J,K)=U(J,K)+V(J,I)*G ENDDOENDDODO J=1,PU(J,K)=U(J,K)/D2ENDDOAPY(K)=APY(K)/D2BY(K)=D2/D1D1=D2ENDDOV(1,1)=1.0V(2,1)=-APY(1)V(2,2)=1.0DO I=1,PDO J=1,QA(I,J)=0.0ENDDOENDDODO I=3,QV(I,I)=V(I-1,I-1)V(I,I-1)=-APY(I-1)*V(I-1,I-1)+V(I-1,I-2)IF(I>=4) THENDO K=I-2,2,-1V(I,K)=-APY(I-1)*V(I-1,K)+V(I-1,K-1)-BY(I-1)*V(I-2,K) ENDDO ENDIFV(I,1)=-APY(I-1)*V(I-1,1)-BY(I-1)*V(I-2,1)ENDDO DO I=1,PIF(I==1) THENT(1)=1.0T1(1)=1.0ELSEIF(I==2) THENT(1)=-APX(1)T(2)=1.0T2(1)=T(1)T2(2)=T(2)ELSET(I)=T2(I-1)T(I-1)=-APX(I-1)*T2(I-1)+T2(I-2) IF(I>=4) THENDO K=I-2,2,-1T(K)=-APX(I-1)*T2(K)+T2(K-1)-BX(I-1)*T1(K) ENDDOENDIFT(1)=-APX(I-1)*T2(1)-BX(I-1)*T1(1)T2(I)=T(I)DO K=I-1,1,-1T1(K)=T2(K)T2(K)=T(K)ENDDOENDIFDO J=1,QDO K=I,1,-1DO L=J,1,-1A(K,L)=A(K,L)+U(I,J)*T(K)*V(J,L) ENDDOENDDOENDDOENDDODT1=0.0DO I=1,NX1=X(I)DO J=1,MY1=Y(J)X2=1.0DD=0.0DO K=1,PG=A(K,Q)DO KK=Q-1,1,-1G=G*Y1+A(K,KK)ENDDOG=G*X2DD=DD+GX2=X2*X1ENDDODT=DD-Z(I,J)DT1=DT1+DT*DTENDDOENDDORETURNEND三、计算结果数表(x,y,f(x,y)): X Y UX TY F(X,Y) 0.00 0.500 1.345 0.243 0.17E+000.00 0.550 1.322 0.269 0.66E+000.00 0.600 1.299 0.295 0.35E+000.00 0.650 1.277 0.322 0.94E+000.00 0.700 1.255 0.350 0.30E-020.00 0.750 1.235 0.377 -0.87E-010.00 0.800 1.215 0.406 -0.58E+000.00 0.850 1.196 0.434 -0.72E+000.00 0.900 1.177 0.463 -0.54E+000.00 0.950 1.159 0.492 -0.86E+000.00 1.050 1.125 0.550 -0.74E+00 0.00 1.100 1.109 0.580 -0.06E+00 0.00 1.150 1.093 0.609 -0.00E+00 0.00 1.200 1.0790.639 -0.18E+00 0.00 1.250 1.064 0.669 -0.52E+00 0.00 1.3001.050 0.699 -0.19E+00 0.00 1.350 1.037 0.729 -0.48E+00 0.001.400 1.024 0.759 -0.68E+00 0.00 1.450 1.011 0.790 -0.52E+00 0.00 1.500 1.000 0.820 -0.29E+000.08 0.500 1.415 0.228 0.67E+00 0.08 0.550 1.391 0.253 0.08E+00 0.08 0.600 1.368 0.279 0.02E+00 0.08 0.650 1.346 0.306 0.47E+00 0.08 0.700 1.325 0.333 0.57E+00 0.08 0.750 1.304 0.360 0.48E-01 0.08 0.800 1.284 0.388 -0.73E-01 0.08 0.850 1.265 0.416 -0.16E+00 0.08 0.900 1.246 0.444 -0.29E+00 0.08 0.950 1.229 0.473 -0.36E+00 0.08 1.000 1.211 0.502 -0.08E+00 0.08 1.050 1.194 0.531 -0.29E+00 0.08 1.100 1.178 0.560 -0.78E+00 0.08 1.150 1.163 0.589 -0.93E+00 0.08 1.200 1.148 0.619 -0.44E+00 0.08 1.250 1.133 0.649 -0.92E+00 0.08 1.300 1.119 0.679 -0.71E+000.08 1.400 1.093 0.739 -0.37E+00 0.08 1.450 1.080 0.769-0.83E+00 0.08 1.500 1.068 0.799 -0.92E+000.16 0.500 1.483 0.214 0.31E+00 0.16 0.550 1.460 0.239 0.64E+00 0.16 0.600 1.437 0.264 0.91E+00 0.16 0.650 1.414 0.290 0.06E+00 0.16 0.700 1.393 0.316 0.70E+00 0.16 0.750 1.372 0.343 0.59E+00 0.16 0.800 1.352 0.370 0.12E+00 0.16 0.850 1.333 0.398 0.77E-02 0.16 0.900 1.315 0.426 -0.83E-01 0.16 0.950 1.297 0.454-0.58E+00 0.16 1.000 1.279 0.483 -0.20E+00 0.16 1.050 1.2620.512 -0.11E+00 0.16 1.100 1.246 0.541 -0.74E+00 0.16 1.1501.231 0.570 -0.09E+00 0.16 1.200 1.216 0.600 -0.59E+00 0.16 1.250 1.201 0.629 -0.66E+00 0.16 1.300 1.187 0.659 -0.71E+00 0.16 1.350 1.174 0.689 -0.32E+00 0.16 1.400 1.161 0.718-0.56E+00 0.16 1.450 1.148 0.748 -0.31E+00 0.16 1.500 1.136 0.778 -0.75E+000.24 0.500 1.551 0.201 0.66E+01 0.24 0.550 1.527 0.2250.03E+000.24 0.650 1.482 0.275 0.64E+00 0.24 0.700 1.460 0.3010.47E+00 0.24 0.750 1.439 0.327 0.34E+00 0.24 0.800 1.419 0.354 0.24E+00 0.24 0.850 1.400 0.381 0.69E+00 0.24 0.900 1.381 0.409 0.04E-01 0.24 0.950 1.363 0.437 -0.42E-01 0.24 1.000 1.346 0.465 -0.06E+00 0.24 1.050 1.329 0.494 -0.59E+00 0.24 1.100 1.313 0.523 -0.83E+00 0.24 1.150 1.297 0.552 -0.15E+00 0.24 1.200 1.282 0.581 -0.19E+00 0.24 1.250 1.267 0.610 -0.84E+00 0.24 1.300 1.253 0.640 -0.66E+00 0.24 1.350 1.240 0.669 -0.30E+00 0.24 1.400 1.227 0.699 -0.86E+00 0.24 1.450 1.214 0.729 -0.84E+00 0.24 1.500 1.202 0.759 -0.77E+000.32 0.500 1.617 0.188 0.28E+01 0.32 0.550 1.593 0.212 0.49E+01 0.32 0.600 1.570 0.236 0.68E+00 0.32 0.650 1.547 0.261 0.75E+00 0.32 0.700 1.526 0.286 0.60E+00 0.32 0.750 1.505 0.312 0.77E+00 0.32 0.800 1.485 0.339 0.05E+00 0.32 0.850 1.466 0.365 0.99E+00 0.32 0.900 1.447 0.393 0.27E+00 0.32 1.000 1.411 0.448 -0.01E-02 0.32 1.050 1.395 0.477-0.41E-01 0.32 1.100 1.378 0.505 -0.18E+00 0.32 1.150 1.3630.534 -0.25E+00 0.32 1.200 1.347 0.563 -0.29E+00 0.32 1.2501.333 0.592 -0.90E+00 0.32 1.300 1.319 0.621 -0.00E+00 0.32 1.350 1.305 0.650 -0.40E+00 0.32 1.400 1.292 0.680 -0.54E+00 0.32 1.450 1.279 0.710 -0.79E+00 0.32 1.500 1.267 0.739-0.91E+000.40 0.500 1.681 0.177 0.91E+01 0.40 0.550 1.658 0.1990.00E+01 0.40 0.600 1.634 0.223 0.83E+01 0.40 0.650 1.612 0.247 0.02E+01 0.40 0.700 1.591 0.272 0.94E+00 0.40 0.750 1.570 0.298 0.49E+00 0.40 0.800 1.550 0.324 0.94E+00 0.40 0.850 1.530 0.350 0.40E+00 0.40 0.900 1.512 0.377 0.33E+00 0.40 0.950 1.493 0.405 0.99E+00 0.40 1.000 1.476 0.432 0.68E+00 0.40 1.050 1.459 0.460 0.08E-01 0.40 1.100 1.443 0.488 -0.84E-01 0.40 1.150 1.427 0.517-0.98E+00 0.40 1.200 1.412 0.545 -0.27E+00 0.40 1.250 1.397 0.574 -0.06E+000.40 1.350 1.369 0.632 -0.66E+00 0.40 1.400 1.356 0.662-0.37E+00 0.40 1.450 1.343 0.691 -0.43E+00 0.40 1.500 1.331 0.721 -0.12E+000.48 0.500 1.745 0.166 0.69E+01 0.48 0.550 1.721 0.188 0.02E+01 0.48 0.600 1.698 0.211 0.74E+01 0.48 0.650 1.676 0.235 0.40E+01 0.48 0.700 1.654 0.259 0.23E+01 0.48 0.750 1.634 0.284 0.56E+00 0.48 0.800 1.613 0.310 0.28E+00 0.48 0.850 1.594 0.336 0.49E+00 0.48 0.900 1.575 0.363 0.31E+00 0.48 0.950 1.557 0.390 0.66E+00 0.48 1.000 1.539 0.417 0.30E+00 0.48 1.050 1.522 0.444 0.34E+00 0.48 1.100 1.506 0.472 0.07E-01 0.48 1.150 1.490 0.500 -0.62E-01 0.48 1.200 1.475 0.529 -0.45E+00 0.48 1.250 1.460 0.557 -0.86E+00 0.48 1.300 1.446 0.586 -0.39E+00 0.48 1.350 1.432 0.615 -0.22E+00 0.48 1.400 1.419 0.644 -0.67E+00 0.48 1.450 1.406 0.674-0.55E+00 0.48 1.500 1.394 0.703 -0.14E+000.56 0.500 1.808 0.156 0.48E+010.56 0.600 1.761 0.200 0.10E+01 0.56 0.650 1.739 0.2230.68E+01 0.56 0.700 1.717 0.247 0.94E+01 0.56 0.750 1.696 0.272 0.33E+01 0.56 0.800 1.676 0.297 0.11E+00 0.56 0.850 1.657 0.323 0.63E+00 0.56 0.900 1.638 0.349 0.97E+00 0.56 0.950 1.620 0.375 0.52E+00 0.56 1.000 1.602 0.402 0.56E+00 0.56 1.050 1.585 0.429 0.47E+00 0.56 1.100 1.568 0.457 0.20E+00 0.56 1.150 1.552 0.485 0.13E+00 0.56 1.200 1.537 0.513 0.09E-01 0.56 1.250 1.522 0.541 -0.47E-01 0.56 1.300 1.508 0.570 -0.99E+00 0.56 1.350 1.4940.599 -0.82E+00 0.56 1.400 1.481 0.627 -0.26E+00 0.56 1.4501.468 0.657 -0.71E+00 0.56 1.500 1.455 0.686 -0.98E+000.64 0.500 1.870 0.147 0.74E+01 0.64 0.550 1.846 0.1680.10E+01 0.64 0.600 1.823 0.190 0.54E+01 0.64 0.650 1.801 0.213 0.42E+01 0.64 0.700 1.779 0.236 0.56E+01 0.64 0.750 1.758 0.260 0.03E+01 0.64 0.800 1.738 0.285 0.42E+01 0.64 0.850 1.718 0.310 0.41E+010.64 0.950 1.681 0.362 0.36E+00 0.64 1.000 1.664 0.388 0.18E+00 0.64 1.050 1.646 0.415 0.28E+00 0.64 1.100 1.630 0.443 0.07E+00 0.64 1.150 1.614 0.470 0.66E+00 0.64 1.200 1.598 0.498 0.09E+00 0.64 1.250 1.584 0.526 0.50E-01 0.64 1.300 1.569 0.554 -0.88E-01 0.64 1.350 1.555 0.583 -0.76E+00 0.64 1.400 1.542 0.611 -0.66E+00 0.64 1.450 1.529 0.640 -0.33E+00 0.64 1.500 1.516 0.669 -0.56E+00 0.72 0.500 1.931 0.139 0.94E+01 0.72 0.550 1.907 0.159 0.84E+01 0.72 0.600 1.884 0.181 0.36E+01 0.72 0.650 1.862 0.203 0.40E+01 0.72 0.700 1.840 0.226 0.47E+01 0.72 0.750 1.819 0.249 0.56E+01 0.72 0.800 1.799 0.273 0.19E+01 0.72 0.850 1.779 0.298 0.37E+01 0.72 0.900 1.760 0.323 0.86E+01 0.72 0.950 1.742 0.349 0.76E+00 0.72 1.000 1.724 0.375 0.24E+00 0.72 1.050 1.707 0.402 0.55E+00 0.72 1.100 1.691 0.429 0.97E+00 0.72 1.150 1.675 0.456 0.27E+00 0.72 1.200 1.659 0.484 0.31E+000.72 1.300 1.630 0.539 0.49E+00 0.72 1.350 1.616 0.5680.72E-02 0.72 1.400 1.602 0.596 -0.69E-01 0.72 1.450 1.589 0.625 -0.67E+00 0.72 1.500 1.576 0.653 -0.20E+000.80 0.500 1.992 0.131 0.31E+01 0.80 0.550 1.968 0.1510.44E+01 0.80 0.600 1.945 0.172 0.41E+01 0.80 0.650 1.922 0.193 0.45E+01 0.80 0.700 1.900 0.216 0.00E+01 0.80 0.750 1.879 0.239 0.10E+01 0.80 0.800 1.859 0.263 0.16E+01 0.80 0.850 1.840 0.287 0.52E+01 0.80 0.900 1.821 0.312 0.02E+01 0.80 0.950 1.802 0.337 0.38E+01 0.80 1.000 1.784 0.363 0.89E+01 0.80 1.050 1.767 0.389 0.28E+00 0.80 1.100 1.751 0.416 0.09E+00 0.80 1.150 1.734 0.4430.23E+00 0.80 1.200 1.719 0.470 0.93E+00 0.80 1.250 1.704 0.498 0.15E+00 0.80 1.300 1.689 0.525 0.86E+00 0.80 1.350 1.675 0.553 0.64E+00 0.80 1.400 1.662 0.582 0.74E-01 0.80 1.450 1.649 0.610 -0.37E-01 0.80 1.500 1.636 0.638 -0.81E+00K和σ分别为:0 0.93E+031 0.61E+012 0.92E-023 0.53E-034 0.16E-055 0.77E-07系数矩阵Crs(按行)为:0.00E+01 -0.83E+01 0.56E+00 0.97E+00 -0.03E+00 0.70E-010.91E+01 -0.99E+00 -0.96E+01 0.17E+01 -0.66E+00 0.10E-01 0.77E+00 0.42E+01 -0.10E+00 -0.81E+00 0.81E+00 -0.62E-01-0.25E+00 -0.21E+00 0.97E+00 -0.18E+00 0.49E+00 -0.63E-010.34E+00 -0.56E+00 0.69E-01 0.51E+00 -0.77E-01 0.27E-01-0.94E-01 0.94E+00 -0.58E+00 0.69E-01 -0.50E-01 0.53E-02 数表(x,y,f(x,y),p(x,y)):X Y F(X,Y) P(X,Y)0.100 0.700 0.58E+00 0.05E+000.100 1.100 -0.66E+00 -0.26E+00 0.100 1.300 -0.68E+00-0.31E+00 0.100 1.500 -0.52E+00 -0.49E+000.200 0.700 0.54E+00 0.19E+00 0.200 0.900 -0.63E-01 -0.65E-01 0.200 1.100 -0.90E+00 -0.90E+00 0.200 1.300 -0.84E+00 -0.90E+00 0.200 1.500 -0.03E+00 -0.04E+000.300 0.700 0.82E+00 0.09E+00 0.300 0.900 0.48E+00 0.11E+00 0.300 1.100 -0.63E+00 -0.88E+00 0.300 1.300 -0.72E+00 -0.96E+00 0.300 1.500 -0.34E+00 -0.84E+000.400 0.700 0.79E+00 0.89E+00 0.400 0.900 0.56E+00 0.63E+00 0.400 1.100 -0.83E-01 -0.04E-01 0.400 1.300 -0.72E+00 -0.71E+00 0.400 1.500 -0.85E+00 -0.07E+000.500 0.700 0.56E+01 0.92E+01 0.500 0.900 0.51E+00 0.23E+00 0.500 1.100 0.59E+00 0.27E+00 0.500 1.300 -0.53E+00 -0.11E+00 0.500 1.500 -0.67E+00 -0.33E+000.600 0.900 0.14E+00 0.75E+00 0.600 1.100 0.19E+00 0.32E+00 0.600 1.300 -0.70E-01 -0.82E-01 0.600 1.500 -0.08E+00 -0.75E+00 0.700 0.700 0.89E+01 0.29E+01 0.700 0.900 0.91E+01 0.11E+010.700 1.100 0.60E+00 0.97E+00 0.700 1.300 0.22E-01 0.06E-01 0.7001.500 -0.53E+00 -0.80E+00 0.800 0.700 0.09E+01 0.06E+01 0.800 0.900 0.32E+01 0.50E+01 0.800 1.100 0.03E+00 0.79E+00 0.800 1.300 0.25E+00 0.50E+00 0.800 1.500 -0.14E+00 -0.28E+00。

北航数值分析大作业一

北航数值分析大作业一

北京航空航天大学数值分析大作业一学院名称自动化专业方向控制工程学号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 积分法计算所需节点数多,计算量大。

数值分析大作业课题一

迭代格式的比较1、自定义函数:function y=fun1(x)y=(3*x+1)/x^2;主程序:x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun1(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果:Please input the initial value:x0=1Please input the number of iterations:k=5x=1.000000x=4.000000x=0.812500x=5.207101x=0.613018x=7.554877结论:发散2、自定义函数:function y=fun2(x)y=(x^3-1)/3;主程序:x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun2(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果:Please input the initial value:x0=1Please input the number of iterations:k=8x=1.000000x=0.000000x=-0.333333x=-0.345679x=-0.347102x=-0.347273x=-0.347294x=-0.347296x=-0.347296结论:收敛3、自定义函数:function y=fun3(x)y=(3*x+1)^(1/3);主程序:x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun3(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果:Please input the initial value:x0=1Please input the number of iterations:k=11x=1.000000x=1.587401x=1.792790x=1.854542x=1.872325x=1.877384x=1.878819x=1.879225x=1.879340x=1.879372x=1.879382x=1.879384结论:收敛4、自定义函数:function y=fun4(x)y=1/(x^2-3);主程序:x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun4(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果:Please input the initial value:x0=1Please input the number of iterations:k=9x=1.000000x=-0.500000x=-0.363636x=-0.348703x=-0.347414x=-0.347306x=-0.347297x=-0.347296x=-0.347296x=-0.347296结论:收敛5、自定义函数:function y=fun5(x)y=sqrt(3+1/x);主程序:x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun5(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果:Please input the initial value:x0=1Please input the number of iterations:k=6x=1.000000x=2.000000x=1.870829x=1.880033x=1.879336x=1.879389x=1.879385结论:收敛6、自定义函数:function y=fun6(x)y=x-1/3*((x^3-3*x-1)/(x^2-1));主程序:x0=input('Please input the initial value£ºx0=');k=input('Please input the number of iterations£ºk='); x(1)=x0;i=1;while i<=kx(i+1)=fun6(x(i));i=i+1;xk=x(i);endfprintf(' x=%f\n',x);运行结果:Please input the initial value:x0=3Please input the number of iterations:k=6x=3.000000x=2.291667x=1.965507x=1.884402x=1.879404x=1.879385x=1.879385结论:收敛要求2,输出迭代次数:自定义函数:interation.m%x0为初始值,k为最大迭代次数,e为控制精度function interation(x0,k,e)x(1)=x0;i=1;while i<=kx(i+1)=fun6(x(i));if abs(x(i+1)-x(i))<efprintf(' the number of iterations: n=%f\n',i);breakendi=i+1;endx=x(i);fprintf(' a root of the original equation: x=%8.6f\n',x);取形式6为例,运行结果:>>interation(5,24,0.00001)the number of iterations: n=7.000000a root of the original equation: x=1.879385。

BUAA数值分析大作业三

北京航空航天大学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 相同。

北航数值分析大作业3(学硕)

《数值分析》作业三院系:机械学院学号: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) 结束。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
相关文档
最新文档