数值实验报告.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)