北航数值分析第二次大作业(最小二乘 simpson rk)

合集下载

北航 数值分析第二次大作业(带双步位移的QR方法)

北航 数值分析第二次大作业(带双步位移的QR方法)

一、算法设计方案:按题目要求,本程序运用带双步位移的QR方法求解给定矩阵的特征值,并对每一实特征值,求解其相应的特征向量。

总体思路:1)初始化矩阵首先需要将需要求解的矩阵输入程序。

为了防止矩阵在后面的计算中被破坏保存A[][]。

2)对给定的矩阵进行拟上三角化为了尽量减少计算量,提高程序的运行效率,在对矩阵进行QR分解之前,先进行拟上三角化。

由于矩阵的QR 分解不改变矩阵的结构,所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。

这里用函数void QuasiTriangularization()来实现,函数形参为double型N维方阵double a[][N]。

3)对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。

用子程序void QR_decomposition()来实现,将Q、R设为形参,然后将计算出来的结果传入Q和R,然后求出RQ乘积。

4)对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛,对QR分解引入双步位移,适当选取位移量,可以避免进行复数运算。

为了进一步减少计算量,在每次进行QR分解之前,先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。

若可以,则得到特征值,同时对矩阵进行降阶处理;若不可以,则进行QR分解。

这里用函数intTwoStepDisplacement_QR()来实现。

这是用来存储计算得到的特征值的二维数组。

考虑到特征值可能为复数,因此将所有特征值均当成复数处理。

此函数中,QR分解部分用子函数void QR_decompositionMk()实现。

这里形参有三个,分别用来传递引入双步位移后的Mk阵,A矩阵,以及当前目标矩阵的维数m。

5)计算特征向量得到特征值后,计算实特征值相应的特征向量。

这里判断所得特征值的虚数部分是否为零。

求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。

因此先初始化矩阵Array,计算(A-λI),再求解方程组。

北航数值分析大作业二(纯原创,高分版)

北航数值分析大作业二(纯原创,高分版)
(R_4 ,I_4 )=( 1.590313458807e+000, 0.000000000000e+000)
(R_5 ,I_5 )=(-1.493147080915e+000, 0.000000000000e+000)
(R_6 ,I_6 )=(-9.891143464723e-001, 1.084758631502e-001)
-0.8945216982
-0.0993313649
-1.0998317589
0.9132565113
-0.6407977009
0.1946733679
-2.3478783624
2.3720579216
1.8279985523
-1.2630152661
0.6790694668
-0.4672150886
6.220134985374e-001
-1.119962139645e-001
-2.521344456568e+000
-1.306189420531e+000
-3.809101150714e+000
8.132800093357e+000
-1.230295627285e+000
-6.753086301215e-001
而其本质就是
1.令 以及最大迭代步数L;
2.若m≤0,则结束计算,已求出A的全部特征值,判断 或 或m≤2是否成立,成立则转3,否则转4;
3.若 ,则得一个特征值 ,m=m-1,降阶;若 ,则计算矩阵:
的特征值得矩阵A的两个特征值,m=m-2,降阶,转2.;
4.若k≤L,成立则令
k=k+1,转2,否则结束计算,为计算出矩阵A的全部特征值;

高中数学北师大版必修3最小二乘估计课后作业Word版含答案

高中数学北师大版必修3最小二乘估计课后作业Word版含答案

§8最小二乘估计一、非标准1.已知回归直线的斜率的估计值是1.23,样本点的中心是(4,5),则线性回归方程是( )A.y=4+1.23xB.y=5+1.23xC.y=0.08+1.23xD.y=1.23+0.08x解析:由已知得b=1.23,=4,=5,于是a=-b=5-1.23×4=0.08,因此线性回归方程为y=1.23x+0.08.答案:C2.用回归直线方程的系数a,b的最小二乘法估计a,b,使函数Q(a,b)的值最小,则Q函数是指( )A.(y i-a-bx i)2B.(y i-a-bx i)C.y i-a-bx iD.(y i-a-bx i)2答案:A3.某地区调查了2~9岁的儿童的身高,由此建立的身高y(cm)与年龄x(岁)的回归模型为y=8.25x+60.13,下列叙述正确的是( )A.该地区一个10岁儿童的身高为142.63cmB.该地区2~9岁的儿童每年身高约增加8.25cmC.该地区9岁儿童的平均身高是134.38cmD.利用这个模型可以准确地预算该地区每个2~9岁儿童的身高解析:由y=8.25x+60.13知斜率的估计值为8.25,说明每增加一个单位年龄,约增加8.25个单位身高,故选B.答案:B4.某产品的广告费用x与销售额y的统计数据如下表:根据上表可得回归方程y=bx+a中的b为9.4,据此模型预报广告费用为6万元时的销售额为( )A.63.6万元B.65.5万元C.67.7万元D.72.0万元解析:∵=3.5,=42,又y=bx+a必过(),∴42=3.5×9.4+a,∴a=9.1.∴线性回归方程为y=9.4x+9.1.∴当x=6时,y=9.4×6+9.1=65.5(万元).答案:B5.已知x与y之间的几组数据如下表:假设根据上表数据所得线性回归方程为y=bx+a,若某同学根据上表中的前两组数据(1,0)和(2,2)求得的直线方程为y=b'x+a',则以下结论正确的是( )A.b>b',a>a'B.b>b',a<a'C.b<b',a>a'D.b<b',a<a'解析:b'==2,a'=0-2×1=-2.x i y i=0+4+3+12+15+24=58,=3.5,=1+4+9+16+25+36=91, ∴b=,a=×3.5=-,∴b<b',a>a'.故选C.答案:C6.下表是某厂1到4用水量y与月份x之间具有线性相关关系,其线性回归方程为y=-0.7x+a,则a的值为.解析:由已知得=2.5,=3.5,因此3.5=-0.7×2.5+a,解得a=5.25.答案:5.257.在一项关于16艘轮船的研究中,船的吨位区间从192t~3246t,船员的人数从5到32,由船员人数y关于吨位x的回归分析得到:y=9.5+0.0062x,假定两艘轮船的吨位相差1000 t,船员平均人数相差,对于最小的船估计的船员人数是,对于最大的船估计的船员人数是.解析:由线性回归方程知船的吨位每增加1000t,则人数增加0.0062×1000≥6(人),又分别令x=192和3246,即可估算船员人数.答案:6 10 298.2014年6月22日,某市物价部门对本市的5家商场的一天销售量及其价格进行调查,5家由数据对应的散点图可知,销售量y与价格x之间有较强的线性相关关系,其线性回归方程是=-3.2x+40,且m+n=20,则其中的n=.解析:(9+9.5+m+10.5+11)=8+(11+n+8+6+5)=6+,线性回归方程一定经过样本中心(),所以6+=-3.2×+40,即3.2m+n=42,由解得故n=10.答案:109.某连锁经营公司所属5个零售店某月的销售额和利润额如下表:万元(1)画出散点图,观察散点图,说明两个变量是否线性相关;(2)用最小二乘法计算利润额y对销售额x的线性回归方程;(3)当销售额为4千万元时,估计利润额的大小.解:(1)散点图如下图.由散点图可以看出变量x,y线性相关.(2)设线性回归方程是y=bx+a,=3.4,=6,所以b==0.5,a=-b=3.4-6×0.5=0.4,即利润额y对销售额x的线性回归方程为y=0.5x+0.4.(3)当销售额为4千万元时,利润额为y=0.5×4+0.4=2.4(百万元).10.(1)利用所给数据求年需求量与年份之间的线性回归方程y=bx+a;(2)利用(1)中所求出的直线方程预测该地2015年的粮食需求量.解:(1)由所给数据看出,年需求量与年份之间是近似直线上升的,下面来求线性回归方程,为此对数据预处理如下:对预处理后的数据,容易算得=0,=3.2.b===6.5,a=-b =3.2.由上述计算结果,知所求线性回归方程为y-257=b(x-2010)+a=6.5(x-2010)+3.2,即y=6.5(x-2010)+260.2.①(2)利用直线方程①,可预测2015年的粮食需求量为6.5×(2015-2010)+260.2=6.5×5+260.2=292.7(万吨).。

数值分析上机实验最小二乘法

数值分析上机实验最小二乘法

---------------------------------------------------------------最新资料推荐------------------------------------------------------数值分析上机实验最小二乘法数值分析实验报告五最小二乘法一、数值分析实验报告五最小二乘法一、题目设有如下数据题目设有如下数据 xj -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 ( )jf x -1.76 0.42 1.2 1.34 1.43 2.25 4.38 -1.76 0.42 1.2 1.34 1.43 2.25 4.38 用三次多项式拟合这组数据,并绘出图形。

二、用三次多项式拟合这组数据,并绘出图形。

二、方法最小二t(f,[xx(1),xx(n)]) 四、结果 save and run 之后:请输入插值节点 as [x1,x2...] [-3 -2 -1 0 1 2 3] 请输入插值节点处对应的函数值 as [f1,f2...] [-1.76 0.42 1.2 1.34 1.43 2.25 4.38] 请输入要求的插值次数m =3 f = 133/100+121469856021/35184372088832*x-8042142191733/450359结果 save and run 之后:请输入插值节点 as [x1,x2...] [-3 -2 -1 0 1 2 3] 请输入插值节点处对应的函数值 as [f1,f2...] [-1.76 0.42 1.2 1.34 1.43 2.25 4.38] 请输入要求的插值次数m =3 f = 133/100+121469856021/35184372088832*x-8042142191733/4503599 627370496*x+1020815915537309/9007199254740992*x9627370496*x+1020815915537309/9007199254740992*x五、拓展:1 / 2最小二乘法计算方法比较简单,是实际中常用的一种方法,但是必须经计算机来实现,如果要保证精度则需要对大量数据进行拟合,计算量很大。

数值分析课件第七章最小二乘逼近共26页文档

数值分析课件第七章最小二乘逼近共26页文档
数值分析课件第七章最小二乘逼近
51、山气日夕佳,飞鸟相与还。 52、木欣欣以向荣, Nhomakorabea涓涓而始流。
53、富贵非吾愿,帝乡不可期。 54、雄发指危冠,猛气冲长缨。 55、土地平旷,屋舍俨然,有良田美 池桑竹 之属, 阡陌交 通,鸡 犬相闻 。
21、要知道对好事的称颂过于夸大,也会招来人们的反感轻蔑和嫉妒。——培根 22、业精于勤,荒于嬉;行成于思,毁于随。——韩愈
23、一切节省,归根到底都归结为时间的节省。——马克思 24、意志命运往往背道而驰,决心到最后会全部推倒。——莎士比亚
25、学习是劳动,是充满思想的劳动。——乌申斯基
谢谢!

北航数值分析第二次大作业--QR分解

北航数值分析第二次大作业--QR分解

《数值分析A》计算实习题目二姓名学号联系方式班级指导教师2012年10月一、算法设计方案整个程序主要分为四个函数,主函数,拟上三角化函数,QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。

因为在最后一个函数中也存在QR分解,所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法,而是在该函数中调用了QR分解函数,这样增强了代码的复用性,减少了程序长度;但由于时间关系,对阵中方法的运算速度没有进行深入研究。

1.为了减少QR分解法应用时的迭代次数,首先对给定矩阵进行拟上三角化处理。

2.对经过拟上三角化处理的矩阵进行QR分解。

3.注意到计算特征值与特征向量的过程首先要应用前面两个函数,于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。

计算过程中,没有采用goto语句,而是根据流程图采用其他循环方式完成了设计,通过对迭代过程的合并,简化了程序的循环次数,最后在计算特征向量的时候采用了列主元高斯消去法。

二、源程序代码#include<stdio.h>#include<math.h>#include<string.h>int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s;double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10]; double X[10][10],Y[10][10],Qt[10][10],M[10][10];double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0}; double epsilon=1e-12;void main(){void Quasiuppertriangular(double A[][10]);void QRdecomposition(double A[][10]);void DoublestepsQR(double A[][10]);int i,j;for(i=0;i<10;i++){for(j=0;j<10;j++){A[i][j]=sin(0.5*(i+1)+0.2*(j+1));Q[i][j]=0;AA[i][j]=A[i][j];}A[i][i]=1.5*cos(2.2*(i+1));AA[i][i]=A[i][i];}Quasiuppertriangular(A); //调用拟上三角化函数printf( "\n A经过拟上三角化矩阵为:\n\n");for(i=0;i<10;i++) //输出拟上三角化矩阵{for(j=0;j<10;j++){printf("%.12e ",A[i][j]); //输出拟上三角化矩阵}printf( "\n\n");}QRdecomposition(A); //调用QR分解函数printf( " 进行QR分解后,R矩阵为:\n\n"); //输出R矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",R[i][j]);}printf( "\n\n");}printf( " Q矩阵为:\n\n"); //输出Q矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",Q[i][j]);}printf( "\n\n");}printf( " RQ矩阵为:\n\n"); //输出RQ矩阵for(i=0;i<10;i++){for(j=0;j<10;j++){printf("%.12e ",RQ[i][j]);}printf( "\n\n");}DoublestepsQR(A); //调用双步位移函数printf( "\n\n 特征值实部依次为:\n\n"); //输出特征值实部for(j=0;j<10;j++){printf("%.12e ",Re[j]);}printf("\n\n 特征值虚部依次为:\n\n "); //输出特征值虚部for(j=0;j<10;j++){printf("%.12e ",Im[j]);}//按行输出特征向量printf( "\n\n 按行输出实特征根相应特征向量为:\n\n");for(i=0;i<10;i++){if(i==1||i==2||i==5||i==6){continue;}for(j=0;j<10;j++){printf("%.12e ",X[i][j]);}printf( "\n\n");}getchar();}//拟上三角化函数void Quasiuppertriangular(double A[][10]) {for(j=0;j<8;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;T[i]=0;W[i]=0;}m=0;for(i=j+2;i<10;i++){if(A[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j+1;i<10;i++){d=d+pow(A[i][j],2);}d=sqrt(d);c=-d;if(A[j+1][j]<=0){c=d;}h=c*(c-A[j+1][j]);U[j+1]=A[j+1][j]-c;for(i=j+2;i<10;i++){U[i]=A[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*A[k][i];}P[i]=P[i]/h;}t=0;for(i=0;i<10;i++){for(k=0;k<10;k++){T[i]=T[i]+U[k]*A[i][k];}T[i]=T[i]/h;t=t+P[i]*U[i];}t=t/h;for(i=0;i<10;i++){W[i]=T[i]-t*U[i];for(k=0;k<10;k++){A[i][k]=A[i][k]-W[i]*U[k]-U[i]*P[k];if(abs(A[i][k])<1e-12){A[i][k]=0;}}}}}//QR分解函数void QRdecomposition(double A[][10]) {for(i=0;i<10;i++){for(j=0;j<10;j++){RQ[i][j]=0;Q[i][j]=0;R[i][j]=A[i][j];}Q[i][i]=1;}for(j=0;j<9;j++){for(i=0;i<10;i++){U[i]=0;P[i]=0;W[i]=0;}m=0;for(i=j+1;i<10;i++){if(R[i][j]!=0){m=m+1;}}if(m==0){continue;}d=0;for(i=j;i<10;i++){d=d+pow(R[i][j],2);}d=sqrt(d);c=-d;if(R[j][j]<=0){c=d;}h=c*(c-R[j][j]);U[j]=R[j][j]-c;for(i=j+1;i<10;i++){U[i]=R[i][j];}for(i=0;i<10;i++){for(k=0;k<10;k++){W[i]=W[i]+U[k]*Q[i][k];}}for(i=0;i<10;i++){for(k=0;k<10;k++){Q[i][k]=Q[i][k]-((W[i]*U[k])/h);}}for(i=0;i<10;i++){for(k=0;k<10;k++){P[i]=P[i]+U[k]*R[k][i];}P[i]=P[i]/h;}for(i=0;i<10;i++){for(k=0;k<10;k++){R[i][k]=R[i][k]-U[i]*P[k];if(abs(R[i][k])<epsilon){R[i][k]=0;}}}}for(i=0;i<10;i++) //计算A(n+1)=RQ {for(j=0;j<10;j++){for(k=0;k<10;k++){RQ[i][j]=RQ[i][j]+R[i][k]*Q[k][j];}}}}//双步位移法计算特征值特征向量函数void DoublestepsQR(double A[][10]){int L=1000,m=9; //定义最大循环次数for(i=0;i<L;i++){for(;m>-1;){if(abs(A[m][m-1])<=epsilon){Re[m]=A[m][m];m=m-1; //降阶if(m==0) //4{Re[0]=A[0][0];break;}if(m==-1){break;}if(m>1){continue;}}b=-A[m][m]-A[m-1][m-1]; //5c=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];if(m==1) //6{if((b*b-4*c)>=0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-1; //循环出口条件break;}if((m>1)&&(abs(A[m-1][m-2])>epsilon)) //8{if(i==L-1){printf("No results! \n");m=0; //循环出口条件break;}break;}if((m>1)&&(abs(A[m-1][m-2])<=epsilon)) //7 {if((b*b-4*c)>0){Re[m]=(-b+sqrt(b*b-4*c))/2;Re[m-1]=(-b-sqrt(b*b-4*c))/2;}if((b*b-4*c)<0){Re[m]=-b/2; Im[m]=sqrt(4*c-b*b)/2;Re[m-1]=-b/2; Im[m-1]=-sqrt(4*c-b*b)/2;}m=m-2; //降阶if(m>0){continue;}if(m==0){Re[0]=A[0][0];break;}}}if(m<=0){break;}s=A[m-1][m-1]+A[m][m]; //9t=A[m][m]*A[m-1][m-1]-A[m][m-1]*A[m-1][m];for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=0;Q[j][k]=0;M[j][k]=0;X[j][k]=0;Y[j][k]=0;}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){M[j][k]=M[j][k]+A[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){M[j][k]=M[j][k]-s*A[j][k];}M[j][j]=M[j][j]+t;}//调用QR分解函数对M矩阵进行分解并传递参数矩阵QQRdecomposition(M);for(j=0;j<10;j++){for(k=0;k<10;k++){Qt[j][k]=Q[k][j];}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){X[j][k]=X[j][k]+Qt[j][l]*A[l][k];}}}for(j=0;j<m+1;j++){for(k=0;k<m+1;k++){for(l=0;l<m+1;l++){Y[j][k]=Y[j][k]+X[j][l]*Q[l][k];}}}for(j=0;j<10;j++){{A[j][k]=Y[j][k];}}}//应用列主元高斯消元法计算实部特征向量for(l=0;l<10;l++){if(l==1||l==2||l==5||l==6){continue;}for(k=0;k<10;k++){for(m=0;m<10;m++){A[k][m]=AA[k][m];}A[k][k]=A[k][k]-Re[l];}for(j=0;j<9;j++){m=j;for(i=j+1;i<10;i++){if(abs(A[i][j])>abs(A[m][j])){m=i;}}{Y[j][k]=A[j][k];A[j][k]=A[m][k];A[m][k]=Y[j][k];}for(k=j+1;k<10;k++){b=A[k][j]/A[j][j];for(i=j;i<10;i++){A[k][i]=A[k][i]-A[j][i]*b;}}}X[l][9]=1;for(i=8;i>=0;i--){c=0;for(j=i+1;j<10;j++){c=c+A[i][j]*X[l][j];}X[l][i]=-c/A[i][i];}}}三、程序输出结果1819。

最小二乘平差 自由度

最小二乘平差 自由度

最小二乘平差最小二乘平差(Least Squares Adjustment)是测量数据处理中的一种常用方法,用于对测量观测数据进行最优估计。

通过该方法可以减小测量误差,提高测量精度,对于工程测量、测绘、地理信息系统等领域具有重要的应用价值。

原理介绍最小二乘平差的原理基于最小化观测量与观测值之间的残差平方和。

在测量过程中,常常会存在观测误差、系统误差等不确定因素,这些因素会导致观测值与真实值之间存在差异。

最小二乘平差通过对所有观测值进行加权,使得观测值与真实值之间的差异最小化。

设有n个观测值,每个观测值的观测量为O,真实值为T,观测误差为e。

则最小二乘平差的目标是找到最优的拟合/估计值x,使得:formula1通过对以上目标进行求解,可以得到最优的拟合/估计值x。

其中,观测值和真实值之间的关系可以通过各种数学模型进行描述,例如线性模型、非线性模型等。

应用场景最小二乘平差在测量数据处理中有广泛的应用。

以下是一些常见的应用场景:1.三角测量:在工程测量中,常用三角测量方法测量不同点之间的距离、角度等。

利用最小二乘平差,可以修正测量误差,提高测量精度。

2.高程测量:在测绘、地理信息系统中,通常需要测量地点的高程信息。

最小二乘平差可以对高程测量数据进行优化处理,提高高程测量的精度。

3.GPS定位:全球卫星定位系统(GPS)在导航、地图绘制等领域有着广泛的应用。

最小二乘平差可用于对GPS观测数据进行处理,提高定位的准确性。

4.建筑变形监测:在建筑工程中,对于建筑物的变形监测和建筑物的结构健康状况评估,最小二乘平差可用于对监测数据进行处理,及时发现异常情况。

实现方法最小二乘平差的实现方法有多种,常用的包括:1.高斯-马尔可夫模型:基于线性模型,通过最小二乘法对观测数据进行拟合估计。

该方法适用于满足高斯分布假设的情况。

2.递归最小二乘法:将观测数据分为多个子集,通过递归的方式对子集数据进行最小二乘拟合,然后合并得到最终的拟合结果。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
3.442724
50248449
2.3
3.787994
64875364
3.787995
11610470
3.787995
11615298
3.787995
11615314
3.787995
11615500
2.4
4.160165
86770901
4.160166
37838632
4.160166
37843906
68954997
2.081044
68957296
2.081044
68957297
2.081044
68957283
1.8
2.308420
90889291
2.308421
16958140
2.308421
16960836
2.308421
16960837
2.308421
16960820
1.9
2.557186
0.3818300505
0.4179591837
0.1012285363
0.2223810345
0.3137066459
0.3626837834
0.0812743884
0.1806481607
0.2606106964
0.3123470770
0.3302393550
截断误差:

最高精度2n-1阶
二、计算结果
2.1
3.123238
20863730
3.123238
59095791
3.123238
59099743
3.123238
59099745
3.123238
59099784
2.2
3.44272
407791075
3.442724
50243946
3.442724
50248333
3.442724
50248341
{
b[i]=0;
for(k=0;k<N;k++)
b[i]+=Base(x,i,k)*y[k];
}
Guass(a,b,c);
s=0;
for(k=0;k<N;k++)
{
yy[k]=0;
for(i=0;i<3;i++)
yy[k]+=c[i]*Base(x,i,k);
s+=pow(yy[k]-y[k],2);
5.25E-04
4点
3.1416119052
1.93E-05
5点
3.1415926966
4.30E-08
6点
3.1415926112
4.24E-08
7点
3.1415926564
2.82E-09
8点
3.1415926538
2.43E-10
程序运行截图:
三、结论
由复化Simpson求积公式的截断误差公式可知其具有4阶收敛性。对于本算例,当步长 时,绝对误差-步长拟合函数关系式近似为:
(2)使用Gauss求积公式重新计算上述积分,比较各阶Gauss积分的数值精度和对 的近似精度。
一、算法描述:
1)Simpson复化求积公式:
截断误差:

截断误差随h减小而减小,然而过小的h可能会造成累积误差。
2)Gauss-Legendre求积公式(权函数 ):
Gauss-Legendre求积节点和求积系数
4.989698
18726518
4.98969
18726587
4.989698
18727057
2.7
5.449311
54277910
5.449312
18554878
5.449312
18561509
5.449312
18561606
5.449312
18562187
2.8
5.940331
96328716
5.940332
}
for(i=0;i<3;i++)
printf("c%d:%f\n",i+1,c[i]);
printf("节点拟合误差:%f\n",s);
}
float Base(float x[],int m,int n)
{
float b;
b=sin((m+1)*x[n]);
return(b);
}
void Guass(float a[][3],float b[],float c[])
1)复化Simpson解:(精确解: )
步长h
1e-1
1e-2
1e-3
1e-4
1e-5
1e-6
1e-7
1e-8
数值积
分解S
3.268
2594
3.154
8593
3.142
9253
3.141
7260
3.141
6060
3.141
5940
3.141
5928
3.141
5926
绝对误
差E
1.27e-1
1.33e-2
60118370
1.686170
60118371
1.686170
60118363
1.6
1.873981
67160519
1.873981
85688338
1.873981
85690253
1.873981
85690254
1.873981
85690243
1.7
2.081044
46724338
2.081044
进而可求解出系数向量 ,从而找到最小二乘拟合函数 。
二、
1)节点数据:
i
1
2
3
4
5
6
7
8
9
9
1
2
3
4
5
6
7
8
9
拟合函数:
精确函数:y=x
拟合误差:251.947067
拟合曲线:
程序运行截图:
2)节点数据:
i
1
2
3
4
5
6
7
8
9
0
0
1.5
3
1.5
0
-1.5
-3
-1.5
0
拟合函数:
拟合误差:0.000000
23964234
2.557186
53989904
2.557186
53993009
2.557186
53993010
2.557186
53992990
2
2.828426
78386010
2.828427
12471086
2.828427
12474610
2.828427
12474611
2.828427
12474587
1.106816
60630836
1.2
1.227879
53964033
1.227879
59356063
1.227879
59356619
1.227879
59356619
1.227879
59356617
1.3
1.364135
90632094ห้องสมุดไป่ตู้
1.364135
99027967
1.364135
99028835
一方面,随着步长的减小,复化Simpson解有效数字增加。另一方面,随着步长的减小,其计算时间逐步增加(步长每缩短10倍,计算时间增长10倍)。因此本算例未考虑步长大于 的情况。理论上,随着 ,应有 。
Gauss-Legendre求积公式收敛速度快,计算精度高。对于本算例,采用5节点Gauss-Legendre求积公式既可使误差控制到 量级,较相同精度的复化Simpson求积法大大减少了计算量。Gauss-Legendre代数精度为2n-1,积分节点数每增加1,代数精度增加2。对于本算例,Gauss-Legendre的有效数字随积分节点数增加而增加,逐步收敛于精确解。
数值分析
第二次上机作业
学号:BY×××
姓名:云××
问题
编写程序编写程序实现如下的线性最小二乘拟合算法:用函数 拟合给定的数据 ,i=0,1,2,…,m。
算法步骤如下:
1.输入m及其数据 ,i=1,2…,m;
2.建立法方程 ;
3.解法方程组;
4.输出 。
用此程序求出如下数据的拟合曲线,计算出拟合误差,并画出这些数据点和拟合曲线的图形,比较分析拟合的效果。
拟合曲线:
程序运行截图:
三、结论
通过比较以上两组给定数据的最小二乘拟合结果可知,拟合基函数对拟合结果的好坏存在很大的影响。对于来自三角函数类型的离散数据点,采用三角函数类型基函数并适当控制基函数个数,可得到良好的拟合结果。
问题
数学上已经证明:
(1)使用Simpson复化求积公式计算 的近似值。选择不同的h,将误差刻画为h的函数绘图,分析该方法的精度。
4.160166
37843935
4.160166
37844207
2.5
4.560358
53227084
4.560359
08668122
4.560359
08673845
4.560359
08673891
4.560359
08674259
2.6
4.989697
相关文档
最新文档