数值实验报告.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;iav+=*（x+i）;
av=av/n;
return（av）;
}
//平方和
floatspfh（intn,float*x）
{inti;
floata,b;
a=0;
for（i=0;ia+=（*（x+i））*（*（x+i））;
return（a）;
}
//和平方
floatshpf（intn,float*x）
{inti;
floata,b;
a=0;
for（i=0;ia=a+*（x+i）;
b=a*a/n;
return（b）;
}
//两数先相乘，再相加
floatdcj（intn,float*x,float*y）
{inti;
floata;
a=0;
for（i=0;ia+=（*（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;ia=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;iprintf（"请输入y的值,请注意与x的值一一对应:
\n"）;
for（i=0;ifor（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;ifor（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;isum+=（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）
升级会员 