数值计算方法编程作业C语言版.docx

数值计算方法编程作业C语言版.docx

《数值计算方法编程作业C语言版.docx》由会员分享，可在线阅读，更多相关《数值计算方法编程作业C语言版.docx（23页珍藏版）》请在冰豆网上搜索。

数值计算方法编程作业C语言版

1:

第二章

（1）二分法求解非线性方程：

#include<>

#include<>

#definef（x）（（x*x-1）*x-1）

voidmain（）

{floata,b,x,eps;

intk=0;

printf（"intputeps\n"）;/*容许误差*/

scanf（"%f",&eps）;

printf（"a,b=\n"）;

for（;;）

{scanf（"%f,%f",&a,&b）;

if（f（a）*f（b）>=0）/*判断是否符合二分法使用的条件*/

printf（"二分法不可使用,请重新输入:

\n"）;

elsebreak;

}

do

{x=（a+b）/2;

k++;

if（f（a）*f（x）<0）/*如果f（a）*f（x）<0，则根在区间的左半部分*/

b=x;

elseif（f（a）*f（x）>0）/*否则根在区间的右半部分*/

a=x;

elsebreak;

}while（fabs（b-a）>eps）;/*判断是否达到精度要求,若没有达到,继续循环*/

x=（a+b）/2;/*取最后的小区间中点作为根的近似值*/

printf（"\nTherootisx=%f,k=%d\n",x,k）;

}

运行结果：

intputeps

a,b=

2,-5

Therootisx=,k=20

Pressanykeytocontinue

总结：

本题关键在于两个端点的取值和误差的判断，此程序较容易。

二分法收敛速度较快，但缺点是只能求解单根。

（2）牛顿法求解非线性方程：

#include<>

#include<>

floatf（floatx）/*定义函数f（x）*/

{return（（-3*x+4）*x-5）*x+6;}

floatf1（floatx）/*定义函数f（x）的导数*/

{return（-9*x+8）*x-5;}

voidmain（）

{floateps,x0,x1=;

printf（"inputeps:

\n"）;

scanf（"%f",&eps）;/*输入容许误差*/

do

{x0=x1;/*准备下一次迭代的初值*/

x1=x0-f（x0）/f1（x0）;/*牛顿迭代*/

}while（fabs（x1-x0）>eps）;/*当满足精度,输出近似根*/

printf（"x=%f\n",x1）;

}

程序运行结果：

x=

总结：

关键是牛顿迭代的应用，程序中最大缺点是函数及其导数已唯一给出确定不可求的随意函数的根，牛顿法比二分法收敛快，可以求重根。

2:

第三章

（1）列主元素消去法求解线性方程：

#include

#include

#defineN20

usingnamespacestd;

voidload（）;

floata[N][N];

intm;

intmain（）{

inti,j;

intc,k,n,p,r;

floatx[N],l[N][N],s,d;

cout<<"下面请输入未知数的个数m=";

cin>>m;

cout<

cout<<"请按顺序输入增广矩阵a:

"<

load（）;

for（i=0;i

{

for（j=i;j

c=（fabs（a[j][i]）>fabs（a[i][i]））?

j:

i;/*找列最大元素*/

for（n=0;n

{s=a[i][n];a[i][n]=a[c][n];a[c][n]=s;}/*将列最大数防在对角线上*/

for（p=0;p

cout<

cout<

for（k=i+1;k

{

l[k][i]=a[k][i]/a[i][i];

for（r=i;r

a[k][r]=a[k][r]-l[k][i]*a[i][r];

}

}

x[m-1]=a[m-1][m]/a[m-1][m-1];

for（i=m-2;i>=0;i--）

{

d=0;

for（j=i+1;j

d=d+a[i][j]*x[j];

x[i]=（a[i][m]-d）/a[i][i];/*求解*/

}

cout<<"该方程组的解为:

"<

for（i=0;i

cout<<"x["<

//system（"pause"）;

return0;

}

voidload（）

{

inti,j;

for（i=0;i

for（j=0;j

cin>>a[i][j];

}

运行结果：

下面请输入未知数的个数m=3

请按顺序输入增广矩阵a:

1234

5108

4692

4692

0

0

该方程组的解为:

x[0]=x[1]=58x[2]=-34Pressanykeytocontinue

总结：

列主元素消去法的目的是为了防止减去一个较小的数时大数淹没小数，而使结果产生较大误差，本程序关键在每次消元时找到相应列中的最大项，然后交换两行位置，在进行计算。

（2）LU分解法求解线性方程：

#include<>

voidsolve（floatl[][100],floatu[][100],floatb[],floatx[],intn）

{inti,j;

floatt,s1,s2;

floaty[100];

for（i=1;i<=n;i++）/*第一次回代过程开始*/

{s1=0;

for（j=1;j

{

t=-l[i][j];

s1=s1+t*y[j];

}

y[i]=（b[i]+s1）/l[i][i];}

for（i=n;i>=1;i--）/*第二次回代过程开始*/

{

s2=0;

for（j=n;j>i;j--）

{

t=-u[i][j];

s2=s2+t*x[j];

}

x[i]=（y[i]+s2）/u[i][i];

}

}

voidmain（）

{floata[100][100],l[100][100],u[100][100],x[100],b[100];

inti,j,n,r,k;

floats1,s2;

for（i=1;i<=99;i++）/*将所有的数组置零，同时将L矩阵的对角值设为1*/

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

{

l[i][j]=0,u[i][j]=0;

if（j==i）l[i][j]=1;

}

printf（"inputn:

\n"）;/*输入方程组的个数*/

scanf（"%d",&n）;

printf（"inputarrayA:

\n"）;/*读取原矩阵A*/

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

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

scanf（"%f",&a[i][j]）;

printf（"inputarrayB:

\n"）;/*读取列矩阵B*/

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

scanf（"%f",&b[i]）;

for（r=1;r<=n;r++）/*求解矩阵L和U*/

{

for（i=r;i<=n;i++）

{

s1=0;

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

s1=s1+l[r][k]*u[k][i];

u[r][i]=a[r][i]-s1;

}

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

{s2=0;

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

s2=s2+l[i][k]*u[k][r];

l[i][r]=（a[i][r]-s2）/u[r][r];

}

}

printf（"arrayL:

\n"）;/*输出矩阵L*/

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

{

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

printf（"%",l[i][j]）;

printf（"\n"）;

}

printf（"arrayU:

\n"）;/*输出矩阵U*/

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

{

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

printf（"%",u[i][j]）;

printf（"\n"）;

}

solve（l,u,b,x,n）;

printf（"解为:

\n"）;

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

printf（"x%d=%f\n",i,x[i]）;

}

运行结果：

inputn:

3

inputarrayA:

223

477

-245

inputarrayB:

31-7

arrayL:

arrayU:

解为:

x1=

x2=

x3=

Pressanykeytocontinue

总结：

关键是把矩阵分解为L、U两个三角矩阵，回代过程比较简单。

3:

第四章

（1）拉格朗日差值多项式;

#include<>

#include<>

#defineMAX100

voidmain（）

{inti,j,k,m,n,N,mi;

floattmp,mx;

floatX[MAX][MAX],Y[MAX],x[MAX],y[MAX],a[MAX];

printf（"\n输入拟合多项式的次数:

\n"）;

scanf（"%d",&m）;

printf（"\n输入给定点的个数n及坐标（x,y）:

\n"）;

scanf（"%d",&N）;

printf（"\n"）;

for（i=0;i

scanf（"%f,%f",&x[i],&y[i]）;

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

{

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

{

tmp=0;

for（k=0;k

tmp=tmp+pow（x[k],（i+j））;

X[i][j]=tmp;

X[j][i]=X[i][j];

}

}

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

{

tmp=0;

for（k=0;k

tmp=tmp+y[k]*pow（x[k],i）;

Y[i]=tmp;

}

for（j=0;j

{

for（i=j+1,mi=j,mx=fabs（X[j][j]）;i<=m;i++）

if（fabs（X[i][j]）>mx）

{

mi=i;

mx=fabs（X[i][j]）;

}

if（j

{

tmp=Y[j];

Y[j]=Y[mi];

Y[mi]=tmp;

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

{

tmp=X[j][k];

X[j][k]=X[mi][k];

X[mi][k]=tmp;

}

}

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

{

tmp=-X[i][j]/X[j][j];

Y[i]+=Y[j]*tmp;

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

X[i][k]+=X[j][k]*tmp;

}

}

a[m]=Y[m]/X[m][m];

for（i=m-1;i>=0;i--）

{

a[i]=Y[i];

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

a[i]-=X[i][j]*a[j];

a[i]/=X[i][i];

}

printf（"\n所求的二次多项式为:

\n"）;

printf（"P（x）=%f",a[0]）;

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

printf（"+（%f）*x^%d",a[i],i）;

}

运行结果：

输入拟合多项式的次数:

5

输入给定点的个数n及坐标（x,y）:

3

1,2

5,3

4,2

所求的二次多项式为:

P（x）=+*x^1+*x^2+*x^3+*x^4+

01934）*x^5Pressanykeytocontinue

总结：

拉格朗日计算公式中，只需要知道各个点即可

4:

第五章

（1）曲线拟合：

#include<>

#include<>

#defineMAX100

voidmain（）

{inti,j,k,m,n,N,mi;

floattmp,mx;

floatX[MAX][MAX],Y[MAX],x[MAX],y[MAX],a[MAX];

printf（"\n输入拟合多项式的次数:

\n"）;

scanf（"%d",&m）;

printf（"\n输入给定点的个数n及坐标（x,y）:

\n"）;

scanf（"%d",&N）;

printf（"\n"）;

for（i=0;i

scanf（"%f,%f",&x[i],&y[i]）;

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

{

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

{

tmp=0;

for（k=0;k

tmp=tmp+pow（x[k],（i+j））;

X[i][j]=tmp;

X[j][i]=X[i][j];

}

}

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

{

tmp=0;

for（k=0;k

tmp=tmp+y[k]*pow（x[k],i）;

Y[i]=tmp;

}

for（j=0;j

{

for（i=j+1,mi=j,mx=fabs（X[j][j]）;i<=m;i++）

if（fabs（X[i][j]）>mx）

{

mi=i;

mx=fabs（X[i][j]）;

}

if（j

{

tmp=Y[j];

Y[j]=Y[mi];

Y[mi]=tmp;

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

{

tmp=X[j][k];

X[j][k]=X[mi][k];

X[mi][k]=tmp;

}

}

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

{

tmp=-X[i][j]/X[j][j];

Y[i]+=Y[j]*tmp;

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

X[i][k]+=X[j][k]*tmp;

}

}

a[m]=Y[m]/X[m][m];

for（i=m-1;i>=0;i--）

{

a[i]=Y[i];

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

a[i]-=X[i][j]*a[j];

a[i]/=X[i][i];

}

printf（"\n所求的二次多项式为:

\n"）;

printf（"P（x）=%f",a[0]）;

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

printf（"+（%f）*x^%d",a[i],i）;

}

输入拟合多项式的次数:

2

输入给定点的个数n及坐标（x,y）:

5

1,2

5,3

2,4

8,3

-1,5

所求的二次多项式为:

P（x）=+*x^1+*x^2Pressanykeytocontinue

5:

第六章

（1）辛普生求积方法：

#include<>

#defineN16/*等分数*/

floatfunc（floatx）

{floaty;

y=（1+x*x）;

return（y）;

}

voidgedianzhi（floaty[],floata,floath）

{inti;

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

y[i]=func（a+i*h）;

}

floatsimpson（floaty[],floath）

{floats,s1,s2;

inti;

s1=y[1];

s2=;

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

{s1+=y[i+1];/*计算奇数项的函数值之和*/

s2+=y[i];/*计算偶数项的函数值之和*/

}

s=y[0]+y[N]+*s1+*s2;

return（s*h/;

}

main（）

{floata,b,h,s,f[N+1];

scanf（"%f,%f",&a,&b）;

h=（b-a）/（float）N;

gedianzhi（f,a,h）;

s=simpson（f,h）;

printf（"s=%f\n",s）;

}

运行结果：

1,3

s=

Pressanykeytocontinue

总结：

辛普生算法是一种积分方法，采用三点法插值，如果h较小的话，误差很小，因为它的插值余项

，辛普生算法比较精确，程序关键是对所取的点的取和，注意

6:

第七章

（1）改进欧拉法求解常微分方程的初值问题

#include<>

floatfunc（floatx,floaty）

{return（y-x）;

}

floateuler（floatx0,floatxn,floaty0,intN）

{floatx,y,yp,yc,h;

inti;

x=x0;

y=y0;

h=（xn-x0）/（float）N;

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

{yp=y+h*func（x,y）;

x=x0+i*h;

yc=y+h*func（x,yp）;

y=（yp+yc）/;

}

return（y）;

}

main（）

{floatx0,xn,y0,e;

intn;

printf（"\ninputn:

\n"）;

scanf（"%d",&n）;

printf（"inputx0,xn:

\n"）;

scanf（"%f,%f",&x0,&xn）;

printf（"inputy0:

\n"）;

scanf（"%f",&y0）;

e=euler（x0,xn,y0,n）;

printf（"y（%f）=%",y0,e）;

}

inputn:

20

inputx0,xn:

1,6

inputy0:

2

y=anykeytocontinue

（2）四阶龙格—库塔法

#include<>

floatfunc（floatx,floaty）

{return（x-y）;

}

floatrunge_kutta（floatx0,floatxn,floaty0,intN）

{floatx,y,y1,y2,h,xh;

floatd1,d2,d3,d4;

inti;

x=x0;

y=y0;

h=（xn-x0）/（float）N;

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

{xh=x+h/2;

d1=func（x,y）;

d2=func（xh,y+h*d1/;

d3=func（xh,y+h*d2/;

d4=func（xh,y+h*d3）;

y=y+h*（d1+2*d2+2*d3+d4）/;

x=x0+i*h;}

return（y）;

}

main（）

{floatx0,xn,y0,e;

intN;

printf（"\ninputn:

\n"）;

scanf（"%d",&N）;

printf（"inputx0,xn:

\n"）;

scanf（"%f,%f",&x0,&xn）;

printf（"inputy0:

\n"）;

scanf（"%f",&y0）;

e=runge_kutta（x0,xn,y0,N）;

printf（"y（%f）=%",y0,e）;

}

inputn:

10

inputx0,xn:

1,2

inputy0:

5

y=anykeytocontinue

2-2Gauss-Seidel方法

#include<>

#include<>

intgsdl（a,b,n,x,eps）

intn;

doublea[],b[],x[],eps;

{inti,j,u,v;

doublep,t,s,q;

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

{u=i*n+i;

p=;

x[i]=;

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

{v=i*n+j;

p=p+fabs（a[v]）;

}

}

if（p>=fabs（a[u]））

{printf（“fail\n”）;

return（-1）;

}

}

p=eps+;

while（p>=eps）

{for（i=0;i<=n-1;i++）

{t=x[i];s=;

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

{if（j!

=i）

{s=s+a[i*n+j]*x[j];

}

x[i]=（b[i]-s）/a[i*n+j];

q=fabs（x[i]-t）/+fabs（x[i]））;

if（q>p）

{p=q;

}

}

}

return

（1）;

}

main（）

{inti;

doubleeps;

staticdoublea[4][4]={{7,2,1,-2}{9,15,3,-2}{-2,-2,11,5}{1,3,2,13}};

staticdoublex[5],b[4]={4,7,-1,0};

eps=;

if（dsdl（a,b,4,x,eps）>0）

{for（i=0;i<=3;i++）

{printf（“x（%d）=%\n”,i,x[i]）;

}

}