数值实验报告.docx
- 文档编号:28648974
- 上传时间:2023-07-19
- 格式:DOCX
- 页数:59
- 大小:332.16KB
数值实验报告.docx
《数值实验报告.docx》由会员分享,可在线阅读,更多相关《数值实验报告.docx(59页珍藏版)》请在冰豆网上搜索。
数值实验报告
实验课程:
数值分析
专业:
信息与计算科学
班级:
10070241
学号:
18
姓名:
付强
中北大学理学院
实验一函数差值方法
【实验目的】
运用数值分析方法及相应的数学软件、函数差值方法解决函数的近似求解。
【实验内容】
对于给定的一元函数
的n+1个节点值
(j=0,1,…n)。
试用拉格朗日公式求其差值多项式或分段二次拉格朗日插值多项式。
数据如下:
x
0.4
0.55
0.65
0.80
0.95
1.05
y
0.41075
0.57815
0.69675
0.90
1.00
1.25382
球五次拉格朗日插值多项式
或分段三次插值多形式,并计算
的值。
(提示:
结果为
)
【实验所使用的仪器设备与软件平台】
电脑、VC++6.0.
【实验方法与步骤】
首先,在理论层面上了结各种插值法的优缺点,理解各种方法的插值原理。
实验一,采取的是拉格朗日插值法,基本公式如下:
其中,
因此,实现此公式的关键步骤在于用C语言循环语句,实现此差值公式,然后运用数组输入具体实验数字,求的所需要的结果。
实验C语言代码如下:
#include
#include
#include
#include
floatlag(intn,float*x,float*y,floatt)
{inti,j;
floatp,s,la;
p=1;la=0;
for(i=0;i {for(j=0;j {if(j! =i) p=p*(t-*(x+j))/(*(x+i)-*(x+j)); } s=*(y+i)*p; p=1; la=la+s; } return(la); } voidmain() {intm,n,k,i; float*x,*y,la,t; printf("请输入给定点的数量n"); scanf("%d",&n); x=(float*)calloc(n,sizeof(float)); if(x==NULL) {printf("内存分配失败"); exit (1); } y=(float*)calloc(n,sizeof(float)); if(y==NULL) {printf("内存分配失败"); exit (1); } x[0]=0.4;x[1]=0.55;x[2]=0.65;x[3]=0.80;x[4]=0.95;x[5]=1.05; y[0]=0.41075;y[1]=0.57815;y[2]=0.69675;y[3]=0.90;y[4]=1.00;y[5]=1.25382; for(i=0;i {printf("x[%d]=%f",i,*(x+i)); printf("y[%d]=%f\n",i,*(y+i)); } for(i=0;;i++) {printf("请输入t的值(退出则输入0)"); scanf("%f",&t); if(t==0) {printf("确认退出请输入0,计算f(0)请按输入任意数字值"); scanf("%d",&m); if(m==0)break; } printf("请输入需要使用的点的个数"); scanf("%d",&k); la=lag(k,x,y,t); printf("插值结果是: \n"); printf("f(%f)=%f\n",t,la); } free(x); free(y); } 【实验结果】 【结果分析与讨论】 实验结果与实验内容中的提示结果完全一样,甚至比实验提示中所给出的精度都要高。 实验二正交多项式与曲线拟合 【实验目的】 用最小二乘法进行曲线拟合,并比较曲线拟合效果,探索拟合曲线的选择与拟合精度间的关系。 【实验内容】 1、编写出legendre、Chebyshev多项式的程序; 2、从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量 与时间 的拟合曲线。 0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 1 1 1 1 1 1 1 1 1 1 1 【实验所使用的仪器设备与软件平台】 电脑、VC++6.0.、matlab 【实验方法与步骤】 一、按照题目要求,首先对数据进行一次拟合,代码如下: #include #include #include #include floataverage(intn,float*x) {inti; floatav; av=0; for(i=0;i av+=*(x+i); av=av/n; return(av); } //平方和 floatspfh(intn,float*x) {inti; floata,b; a=0; for(i=0;i a+=(*(x+i))*(*(x+i)); return(a); } //和平方 floatshpf(intn,float*x) {inti; floata,b; a=0; for(i=0;i a=a+*(x+i); b=a*a/n; return(b); } //两数先相乘,再相加 floatdcj(intn,float*x,float*y) {inti; floata; a=0; for(i=0;i a+=(*(x+i))*(*(y+i)); return(a); } //两数先相加,再相乘 floatdjc(intn,float*x,float*y) {inti; floata=0,b=0; for(i=0;i {a=a+*(x+i); b=b+*(y+i); } a=a*b/n; return(a); } //系数a floatxsa(intn,float*x,float*y) {floata,b,c,d,e; a=spfh(n,x); b=shpf(n,x); c=dcj(n,x,y); d=djc(n,x,y); e=(c-d)/(a-b); //printf("%f%f%f%f",a,b,c,d); return(e); } floathe(intn,float*y) {inti; floata; a=0; for(i=0;i a=a+*(y+i); return(a); } floatxsb(intn,float*x,float*y,floata) {floatb,c,d; b=he(n,y); c=he(n,x); d=(b-a*c)/n; return(d); } intmain() {intn,i; float*x,*y,a,b; printf("请输入将要输入的有效数值组数n的值: "); scanf("%d",&n); x=(float*)calloc(n,sizeof(float)); if(x==NULL) {printf("内存分配失败"); exit (1); } y=(float*)calloc(n,sizeof(float)); if(y==NULL) {printf("内存分配失败"); exit (1); } printf("请输入x的值\n"); for(i=0;i printf("请输入y的值,请注意与x的值一一对应: \n"); for(i=0;i for(i=0;i {printf("x[%d]=%3.2f",i,*(x+i)); printf("y[%d]=%3.2f\n",i,*(y+i)); } a=xsa(n,x,y); b=xsb(n,x,y,a); printf("经最小二乘法拟合得到的一元线性方程为: \n"); printf("f(x)=%3.2fx+%3.2f\n",a,b); } 二、经过matlab进行曲线拟合可知,该组数据近似的服从3次曲线,因此下面对所给数据进行3次最小二乘法拟合,代码如下: #include #include #include #include #defineN12//N个点 #defineT3//T次拟合 #defineW1//权函数 #definePRECISION0.00001 floatpow_n(floata,intn) { inti; if(n==0) return (1); floatres=a; for(i=1;i { res*=a; } return(res); } voidmutiple(floata[][N],floatb[][T+1],floatc[][T+1]) { floatres=0; inti,j,k; for(i=0;i for(j=0;j { res=0; for(k=0;k { res+=a[i][k]*b[k][j]; c[i][j]=res; } } } voidmatrix_trans(floata[][T+1],floatb[][N]) { inti,j; for(i=0;i { for(j=0;j { b[j][i]=a[i][j]; } } } voidinit(floatx_y[][2],intn) { inti; printf("请输入%d个已知点: \n",N); for(i=0;i { printf("(x%dy%d): ",i,i); scanf("%f%f",&x_y[i][0],&x_y[i][1]); } } voidget_A(floatmatrix_A[][T+1],floatx_y[][2],intn) { inti,j; for(i=0;i { for(j=0;j { matrix_A[i][j]=W*pow_n(x_y[i][0],j); } } } voidprint_array(floatarray[][T+1],intn) { inti,j; for(i=0;i { for(j=0;j { printf("%-g",array[i][j]); } printf("\n"); } } voidconvert(floatargu[][T+2],intn) { inti,j,k,p,t; floatrate,temp; for(i=1;i { for(j=i;j { if(argu[i-1][i-1]==0) { for(p=i;p { if(argu[p][i-1]! =0) break; } if(p==n) { printf("方程组无解! \n"); exit(0); } for(t=0;t { temp=argu[i-1][t]; argu[i-1][t]=argu[p][t]; argu[p][t]=temp; } } rate=argu[j][i-1]/argu[i-1][i-1]; for(k=i-1;k { argu[j][k]-=argu[i-1][k]*rate; if(fabs(argu[j][k])<=PRECISION) argu[j][k]=0; } } } } voidcompute(floatargu[][T+2],intn,floatroot[]) { inti,j; floattemp; for(i=n-1;i>=0;i--) { temp=argu[i][n]; for(j=n-1;j>i;j--) { temp-=argu[i][j]*root[j]; } root[i]=temp/argu[i][i]; } } voidget_y(floattrans_A[][N],floatx_y[][2],floaty[],intn) { inti,j; floattemp; for(i=0;i { temp=0; for(j=0;j { temp+=trans_A[i][j]*x_y[j][1]; } y[i]=temp; } } voidcons_formula(floatcoef_A[][T+1],floaty[],floatcoef_form[][T+2]) { inti,j; for(i=0;i { for(j=0;j { if(j==T+1) coef_form[i][j]=y[i]; else coef_form[i][j]=coef_A[i][j]; } } } voidprint_root(floata[],intn) { inti,j; printf("%d个点的%d次拟合的多项式系数为: \n",N,T); for(i=0;i { printf("a[%d]=%g,",i+1,a[i]); } printf("\n"); printf("拟合曲线方程为: \ny(x)=%g",a[0]); for(i=1;i { printf("+%g",a[i]); for(j=0;j { printf("*X"); } } printf("\n"); } voidprocess() { floatx_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_formu[T+1][T+2],y[T+1],a[T+1]; init(x_y,N); get_A(matrix_A,x_y,N); printf("矩阵A为: \n"); print_array(matrix_A,N); matrix_trans(matrix_A,trans_A); mutiple(trans_A,matrix_A,coef_A); printf("法矩阵为: \n"); print_array(coef_A,T+1); get_y(trans_A,x_y,y,T+1); cons_formula(coef_A,y,coef_formu); convert(coef_formu,T+1); compute(coef_formu,T+1,a); print_root(a,T+1); } voidmain() { process(); } 【实验结果】 【结果分析与讨论】 对一次拟合做matlab检验,结果如下: 可从图像得一次拟合是错误的,于是进行三次拟合,结果如下: 因此,得出上述数据应进行三次最小二乘法拟合才能得到比较精确的解。 实验三数值积分与数值微分 【实验目的】 选用复合梯形公式,复合辛普森公式, 算法,高斯算法计算所给出的积分的近似解。 【实验内容】 选用复合梯形公式,复合辛普森公式, 算法,高斯算法计算 (1) (2) (3) 【实验所使用的仪器设备与软件平台】 电脑、VC++6.0。 【实验方法与步骤】 一、运用 算法,对 (1)进行数值法几分计算代码如下: #include #include staticdoubleT[20][20]; doublef(doublex){ returnpow((4-pow(sin(x),2)),0.5); } doublefuhe(doublea,doubleb,intn){ doublesum=0,h=0,k=0; inti; if(n==0){ h=b-a; returnh*(f(a)+f(b))/2.0; } else{ k=pow(2,n); h=(b-a)/k; for(i=0;i sum+=(f(a+(b-a)*i/k)+f(a+(b-a)*(i+1)/k))*h/2.0; } returnsum; } } doubleromberg(doublea,doubleb,doublee){ intj=0,k=0,i=2; T[0][1]=fuhe(a,b,0); T[1][1]=fuhe(a,b,1); T[0][2]=(4*T[1][1]-T[0][1])/3.0; for(;fabs(T[0][i]-T[0][i-1])>=e;i++){ T[i][1]=fuhe(a,b,i+1); j=2; for(k=i-1;k>=0;k--){ T[k][j]=((pow(4,j-1)*T[k+1][j-1]-T[k][j-1]))/(pow(4,j-1)-1); j++; } } returnT[0][i]; } intmain(){ floata,b,e; doubleresult; printf("请输入积分下限a: \n"); scanf("%f",&a); printf("请输入积分上限b: \n"); scanf("%f",&b); printf("输入e的近似值: \n"); scanf("%f",&e); result=romberg(a,b,e); printf("Basedonyourinputparameterstheresultsis: \n%f\n\n",result); scanf("%f",&e); } 二、运用复合辛普森算法,对 (1)进行数值法几分计算代码如下: #include #include doublef(doublex); doublezdyf(doublex); doubleff(doublea,doubleb,doubleeps); doubleff(doublea,doubleb,doubleeps,double(*f)(double)); voidmain() {charstr[1]; doublea,b,eps,t;//[a,b]内,当前精度eps,接受结果到t; double(*fff)(double); a=0.0; b=0.25; eps=0.0000001; t=ff(a,b,eps); printf("被积函数: \n"); printf("pow((4-pow(sin(x),2)),0.5)\n"); printf("积分运算: \n"); printf("区间[%3.2f,%3.2f]\n",a,b); printf("积分结果: \n"); printf("t=%e\n",t); printf("\n"); gets(str); } /*--------------------自定义被积函数-------------*/ doublef(doublex){ doubley; y=pow((4-pow(sin(x),2)),0.5); return(y); } /*--------------------积函数(辛普森)在区间[a,b]积f(x)-------------*/ doubleff(doublea,doubleb,doubleeps){ intn,k; doubleh,t1,t2,s1,s2,ep,p,x; n=1;//当前子节点个数-1,不包括两端a,b h=b-a;//当前区间长度 t1=h*(f(a)+f(b))/2.0;//计算原始梯形面积 s1=t1;//将原始梯形面积->s1 ep=eps+1.0;//误差初始化 /*-----------------辛普森算法------------*/ while(ep>=eps){//循环比较(当前ep)和(标准eps) p=0.0;//累加暂存量初始化 for(k=0;k<=n-1;k++){//循环求和实现 x=a+(k+0.5)*h;//节点定位 p+=f(x);//节点函数值累加 } t2=(t1+h*p)/2.0;//用(上一个t1)构造(当前t2); s2=(4.0*t2-t1)/3.0;//用(上一个t1)和(当前t2)构造(当前s2) ep=fabs(s2-s1);//用(当前s2)和(上一个s1)构造(当前ep) t1=t2;//为构造下一个t准备 s1=s2;//为构造下一个s准备 n+=n;//节点个数倍增 h/=2.0;//区间二分化 } return(s2);//返回(当前s2) } doubleff(doublea,doubleb,doubleeps,double(*f)(double)) { intn,k; doubleh,t1,t2,s1,s2,ep,p,x; n=1;//当前子节点个数-1,不包括两端a,b h=b-a;//当前区间长度 t1=h*((*f)(a)+(*f)(b))/2.0;//计算原始梯形面积 s1=t1;//将原始梯形面积->s1 ep=eps+1.0;//误差初始化 /*-----------------辛普森算法------------*/ while(ep>=eps){//循环比较(当前ep)和(标准eps)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 实验 报告