13年数值分析试题.docx

13年数值分析试题.docx

《13年数值分析试题.docx》由会员分享，可在线阅读，更多相关《13年数值分析试题.docx（25页珍藏版）》请在冰点文库上搜索。

13年数值分析试题

福建农林大学

研究生课程考试卷

授课时间2012—2013学年度第1学期

学号：

1121213004

研究生姓名：

吴翠兰

课程名称：

数值分析

考试时间：

2012-12-25

考生成绩：

授课或主考教师：

（签章）

2012-13学年

（1）福建农林大学机电学院研究生《数值分析》课程试卷

姓名：

吴翠兰班级：

2012级学号：

1121213004成绩：

1、用对分区间法求方程

在区间[1,2]上的一个实根。

精度要求为

。

（10分）

【解】/*头文件包含区*/

#include

#include

#include

/*定义声明区*/

intN=0,n=0;

doublea=0,b=0,e=0,x=0;

voidstep_1（void）;

voidstep_2（void）;

voidstep_3（void）;

voidstep_4（void）;

voidstep_5（void）;

voidstep_6（void）;

voidstep_7（void）;

/*f（x）函数*/

doublef（doubleX）

{

doubleresult;

//输入f（x）

result=X+log（X）-2.2;

return（result）;

}

/*子函数*/

voidstep_1（void）

{

doublef（doubleX）;

a=1;

b=2;

e=0.000001;

N=1000;

}

voidstep_2（void）

{

n=0;

}

voidstep_3（void）

{

x=（a+b）/2;

}

voidstep_4（void）

{

if（fabs（f（x））<0.000001||（（b-a）/2）

\n"）;printf（"方程的根为\nx=%.5f",x）;;while

（1）;}

elsestep_5（）;

}

voidstep_5（void）

{

if（n

else{printf（"求解失败...\n"）;while

（1）;}

}

voidstep_6（void）

{

if（（f（a）*f（x））<0）b=x;

elsea=x;

}

voidstep_7（void）

{

n=n+1;step_3（）;

}

/*主函数*/

voidmain（）

{

step_1（）;

step_2（）;

while

（1）

{

step_3（）;

step_4（）;

step_7（）;

}

}

调试结果：

2、通过对系数矩阵的三角分解求解如下线性方程组。

（10分）

【解】

/*头文件包含区*/

#include

#include

#include

#defineN3//矩阵阶数：

3

/*定义声明区*/

inti=0,j=0,k=0,r=0;

floatsum1=0,sum2=0,sum3=0,sum4=0;

floatA[N+1][N+1]={0};

floatL[N+1][N+1]={0};

floatU[N+1][N+1]={0};

floatb[N+1]={0};

floaty[N+1]={0};

floatx[N+1]={0};

/*主函数*/

voidmain（）

{

/*矩阵A[i][j]*//*矩阵b[k]*/

A[1][1]=2;A[1][2]=-1;A[1][3]=-1;b[1]=4;

A[2][1]=3;A[2][2]=4;A[2][3]=-2;b[2]=11;

A[3][1]=3;A[3][2]=-2;A[3][3]=4;b[3]=11;

//求矩阵U&矩阵L

for（j=1;j<=N;j++）

{

U[1][j]=A[1][j];

}

for（i=2;i<=N;i++）

{

L[i][1]=A[i][1]/U[1][1];

}

for（r=2;r<=N;r++）

{

for（j=r;j<=N;j++）

{

for（k=1;k<=（r-1）;k++）

{

sum1=sum1+L[r][k]*U[k][j];

}

U[r][j]=A[r][j]-sum1;

sum1=0;

}

for（i=r+1;i<=N;i++）

{

for（k=1;k<=（r-1）;k++）

{

sum2=sum2+L[i][k]*U[k][r];

}

L[i][r]=（A[i][r]-sum2）/U[r][r];

sum2=0;

}

}

for（i=1;i<=N;i++）{L[i][i]=1;}

//求矩阵y

y[1]=b[1];

for（k=2;k<=N;k++）

{

for（j=1;j<=（k-1）;j++）

{

sum3=sum3+L[k][j]*y[j];

}

y[k]=b[k]-sum3;

sum3=0;

}

//求矩阵x

x[N]=y[N]/U[N][N];

for（k=N-1;k>=1;k--）

{

for（j=k+1;j<=N;j++）

{

sum4=sum4+U[k][j]*x[j];

}

x[k]=（y[k]-sum4）/U[k][k];

sum4=0;

}

/*****************************/

/*****************************/

//输出矩阵A

printf（"输出矩阵A:

\n"）;

for（i=1;i<=N;i++）

{

for（j=1;j<=N;j++）

{

printf（"%.3f",A[i][j]）;

}

printf（"\n"）;

}

//输出矩阵L

printf（"输出矩阵L:

\n"）;

for（i=1;i<=N;i++）

{

for（j=1;j<=N;j++）

{printf（"%.3f",L[i][j]）;

}

printf（"\n"）;

}

//输出矩阵U

printf（"输出矩阵U:

\n"）;

for（i=1;i<=N;i++）

{

for（j=1;j<=N;j++）

{

printf（"%.3f",U[i][j]）;

}

printf（"\n"）;

}

//输出矩阵y

printf（"输出矩阵y:

\n"）;

for（k=1;k<=N;k++）

{

printf（"y%d=%.3f",k,y[k]）;

printf（"\n"）;

}

//输出矩阵x

printf（"输出矩阵x:

\n"）;

for（k=1;k<=N;k++）

{printf（"x%d=%.3f",k,x[k]）;

printf（"\n"）;

}

}

调试结果：

3、给定函数表如下，用抛物线插值法求f（x）在x=1.682，x=1.813处的近似值。

（10分）

x

1.615

1.634

1.702

1.828

1.921

y=f（x）

2.41450

2.46459

2.65271

3.03035

3.34066

【解】

#include

#include

#defineN5

typedefstruct

{doublex;doubley;}

Point;

doubleparabola（Point*,Point*,Point*,double）;

intmain（）

{doublex=0,y=0;

inti=0,n=0;//n为要插的中间那个点的位置

charc='y';Pointp[N];

for（i=0;i

{printf（"Pleaseinputthedata%d:

",i+1）;

scanf（"%lf%lf",&p[i].x,&p[i].y）;

while（10!

=getchar（））

{continue;

}

}while（'n'!

=c）

{printf（"Pleaseinputx:

"）;

scanf（"%lf",&x）;

while（10!

=getchar（））

{continue;

}//计算插值点的位子

if（2*x<（p[1].x+p[2].x））n=1;

elseif（2*x>（p[N-2].x+p[N-1].x））n=N-2;

elsen=（int）N/2;//printf（"%d\n",n）;

y=parabola（&p[n-1],&p[n],&p[n+1],x）;

printf（"y=%lf\n",y）;

printf（"Continue?

（y/n）"）;

scanf（"%c",&c）;

while（10!

=getchar（））

{continue;

}

}return0;

}

doubleparabola（Point*p1,Point*p2,Point*p3,doublex）

{doubletemp1=0,temp2=0,temp3=0,y=0;

doublea=0,b=0,c=0;

temp1=（p2->y-p1->y）*（x-p1->x）;

temp2=（p3->y-p1->y）*（p2->x-p1->x）-（p2->y-p1->y）*（p3->x-p1->x）;

temp3=（p3->x-p1->x）*（p3->x-p1->x）*（p2->x-p1->x）;

y=p1->y+temp1/（p2->x-p1->x）+temp2*（x-p1->x）*（x-p2->x）/temp3;returny;

}

调试结果：

4、对于函数

，从x0=0开始，去步长h=0.1的20个数据点，求5次最小二乘

拟合多项式

其中

（即求系数ai,i=0,1,…,5）（10分）

【解】

/*头文件包含区*/

#include

#include

#include

#definem5//最小二乘m次拟合多项式：

5

#definen20//数据点的个数：

20

/*定义声明区*/

inti=0,j=0,k=0,r=0;

doublesum1=0,sum2=0,sum3=0,sum4=0;

doubleA[m+2][m+2]={0};

doubleL[m+2][m+2]={0};

doubleU[m+2][m+2]={0};

doubleB[m+2]={0};

doubleY[m+2]={0};

doubleX[m+2]={0};

doublea[m+1]={0};

doubley[n+1]={0};

doublex[n+1]={0};

/*主函数*/

voidmain（）

{

/*数组x[i]&y[i]*/

/**************************************************************************/

for（i=1;i<=20;i++）

{

x[i]=-0.95+0.1*（i-1）;

}

//x[1]=0;x[2]=0.1;x[3]=0.2;x[4]=0.3;x[5]=0.4;x[6]=0.5;x[7]=0.6;x[8]=0.7;x[9]=0.8;x[10]=0.9;x[11]=1.0;x[12]=1.1;x[13]=1.2;x[14]=1.3;x[15]=1.4;x[16]=1.5;x[17]=1.6;x[18]=1.7;x[19]=1.8;x[20]=1.9;

for（i=1;i<=20;i++）

{

y[i]=x[i]+0.95-exp（-x[i]-0.95）;

}

//y[1]=-0.2209;y[2]=0.3295;y[3]=0.8826;y[4]=1.4392;y[5]=2.0003;y[6]=2.5645;y[7]=3.1334;y[8]=3.7601;y[9]=4.2836;

/**************************************************************************/

//求矩阵A

for（i=1;i<=（m+1）;i++）

{

for（j=1;j<=（m+1）;j++）

{

for（k=1;k<=n;k++）

{

A[i][j]=A[i][j]+pow（x[k],i+j-2）;

}

}

}

//求矩阵B

for（i=1;i<=（m+1）;i++）

{

for（k=1;k<=n;k++）

{

B[i]=B[i]+y[k]*pow（x[k],i-1）;

}

}

//求矩阵U&矩阵L

for（j=1;j<=（m+1）;j++）

{

U[1][j]=A[1][j];

}

for（i=2;i<=（m+1）;i++）

{

L[i][1]=A[i][1]/U[1][1];

}

for（r=2;r<=（m+1）;r++）

{

for（j=r;j<=（m+1）;j++）

{

for（k=1;k<=（r-1）;k++）

{

sum1=sum1+L[r][k]*U[k][j];

}

U[r][j]=A[r][j]-sum1;

sum1=0;

}

for（i=r+1;i<=（m+1）;i++）

{

for（k=1;k<=（r-1）;k++）

{

sum2=sum2+L[i][k]*U[k][r];

}

L[i][r]=（A[i][r]-sum2）/U[r][r];

sum2=0;

}

}

for（i=1;i<=（m+1）;i++）{L[i][i]=1;}

//求矩阵Y

Y[1]=B[1];

for（k=2;k<=（m+1）;k++）

{

for（j=1;j<=（k-1）;j++）

{

sum3=sum3+L[k][j]*Y[j];

}

Y[k]=B[k]-sum3;

sum3=0;

}

//求矩阵X

X[（m+1）]=Y[（m+1）]/U[（m+1）][（m+1）];

for（k=（m+1）-1;k>=1;k--）

{

for（j=k+1;j<=（m+1）;j++）

{

sum4=sum4+U[k][j]*X[j];

}

X[k]=（Y[k]-sum4）/U[k][k];

sum4=0;

}

//求矩阵a

for（i=0;i<=m;i++）{a[i]=X[i+1];}

//输出矩阵A

printf（"输出矩阵A:

\n"）;

for（i=1;i<=（m+1）;i++）

{

for（j=1;j<=（m+1）;j++）

{

printf（"%.3f",A[i][j]）;

}

printf（"\n"）;

}

//输出矩阵B

printf（"输出矩阵B:

\n"）;

for（k=1;k<=（m+1）;k++）

{

printf（"B%d=%.3f",k,B[k]）;

printf（"\n"）;

}

//输出矩阵a

printf（"输出矩阵a:

\n"）;

for（k=0;k<=m;k++）

{

printf（"a%d=%.3f",k,a[k]）;

printf（"\n"）;

}

//输出最小二乘m次拟合多项式

printf（"最小二乘%d次拟合多项式为:

\n",m）;

printf（"f（x）="）;

printf（"%.3f+",a[0]）;

printf（"%.3f*x",a[1]）;

for（k=2;k<=m;k++）

{

if（a[k]>=0）{printf（"+"）;}

printf（"%.3f*x^%d",a[k],k）;

printf（"\n\n"）;

}

}

调试结果：

5、用Simpson公式逐次分半算法计算定积分。

要求误差不超过10-6。

（10分）

（所需节点处的函数值按需自行计算）

【解】

#include

#include

doublefsimpf（doublex）

{

return（log（1+x））/（1+x*x）;

}

doublefsimp（doublea,doubleb,doubleeps）

{

intn,k;

doubleh,t1,t2,s1,s2,ep,p,x;

n=1;h=b-a;

t1=h*（fsimpf（a）+fsimpf（b））/2.0;

s1=t1;

ep=eps+1.0;

while（ep>=eps）

{

p=0.0;

for（k=0;k<=n-1;k++）

{

x=a+（k+0.5）*h;

p=p+fsimpf（x）;

}

t2=（t1+h*p）/2.0;

s2=（4.0*t2-t1）/3.0;

ep=fabs（s2-s1）;

t1=t2;s1=s2;n=n+n;h=h/2.0;

}

return（s2）;

}

voidmain（）

{

doublea,b,eps,t;

a=0.0;b=1.0;eps=0.000001;

//adefiniteintegralbySimpsonMethod.

t=fsimp（a,b,eps）;

printf（"%e\n",t）;

printf（"\nPressanykeytoquit..."）;

getch（）;

}

调试结果：

6、用Admas预测-校正公式求初值问题的数值解。

取步长h=0.1,计算结果取8位小数。

（10分）

【解】

/*头文件包含区*/

#include

#include

#include

/*定义声明区*/

intn;

doublea=0,b=0,h=0,x=0,y=0,N=0,K1=0,K2=0,K3=0,K4=0,x0=0,yo=0;

doubleX[11];

doubleYY[11];

doubleyy[11];

/*f（x,y）函数*/

doublef（doublex,doubley）

{

doubleresult;

//输入f（x,y）

result=1-（2*x+y）/（1+pow（x,2））;

return（result）;

}

/*主函数*/

voidmain（）

{

//输入a,b,h,yo

a=0;

b=2;

h=0.1;

yo=0;

N=（b-a）/h;

n=1;

x0=a;

printf（"数值解（x,y）为:

\n"）;

printf（"（%.1f,%.8f）\n",x,y）;

while

（1）

{

x=x0+h;

K1=f（x0,yo）;

K2=f（x0+h/2,yo+h/2*K1）;

K3=f（x0+h/2,yo+h/2*K2）;

K4=f（x0+h,yo+h*K3）;

y=yo+h/6*（K1+2*（K2+K3）+K4）;

//输出（x,y）

printf（"（%.1f,%.8f）\n",x,y）;

if（n<10）

{

n=n+1;

x0=x;

yo=y;

}

else

break;

}

X[0]=0;

X[1]=0.2;

X[2]=0.4;

X[3]=0.6;

X[4]=0.8;

X[5]=1.0;

X[6]=1.2;

X[7]=1.4;

X[8]=1.6;

X[9]=1.8;

X[10]=2.0;

YY[0]=0;

YY[1]=0.14484138;

YY[2]=0.20285825;

YY[3]=0.21031166;

for（n=3;n<=20;n++）

{

yy[n+1]=YY[n]+（h/24）*（55*f（X[n],YY[n]）-59*f（X[n-1],YY[n-1]）+37*f（X[n-2],YY[n-2]）-9*f（X[n-3],YY[n-3]））;

YY[n+1]=YY[n]+（h/24）*（9*f（X[n+1],yy[n+1]）+19*f（X[n],YY[n]）-5*f（X[n-1],YY[n-1]）+f（X[n-2],YY[n-2]））;

}

for（n=1;n<=20;n++）

{

printf（"YY[%.d]=",n）;

printf（"%.8f\n",YY[n]）;

}

}

调试结果：