数值计算C语言常用小程序DOC.docx

数值计算C语言常用小程序DOC.docx

《数值计算C语言常用小程序DOC.docx》由会员分享，可在线阅读，更多相关《数值计算C语言常用小程序DOC.docx（20页珍藏版）》请在冰豆网上搜索。

数值计算C语言常用小程序DOC

1、秦九韶算法2、二分法3、拉格朗日插值4、埃特金算法5、复化梯形法

6、复化辛甫生算法7、二阶龙格库塔方法8、四阶龙格库塔方法　9、改进的欧拉方法　10、迭代法　11、埃特金加速方法：

12、牛顿迭代法　13、追赶法　14、雅克比迭代　15、蛋白质设计：

17高斯消去法:

1、秦九韶算法

P11.3利用秦九韶算法求多项式，在x=3时的值。

程序：

#include

#include

voidmain（）

{floata[100],v,x;

intn,i,k;

scanf（"%d%f",&n,&x）;

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

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

v=a[n];

k=1;

do

{v=x*v+a[n-k];

k=k+1;

}

while（k<=n）;

printf（"v=%f",v）;}

运行结果：

2、二分法

P11.1用二分法求方程法x*x*x-x-1=0在[1,2]内的近似根，要求误差不超过

#include

#include

floatfun（float）;

voidmain（）

{floata,b,c,x,y,y1;

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

y1=fun（a）;

do

{x=（a+b）/2;

y=fun（x）;

{if（y*y1>0）

a=x;

else

b=x;}}

while（（b-a）>=c）;

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

}

floatfun（floatm）

{floatn;

n=m*m*m-m-1;

return（n）;

}

运行结果：

3、拉格朗日插值

程序：

#include

main（）

{floata,b,t,x[100],y[100];

intn,i,j,k;

scanf（"%f%d",&a,&n）;

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

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

k=0;

b=0;

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

{t=1;

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

{if（j!

=k）

t=t*（a-x[j]）/（x[k]-x[j]）;}

b=b+t*y[k];

}

printf（"%f\n",b）;

}

4、埃特金算法

程序：

#include

#include

main（）

{floata,b,c,x[100],y[100];

inti,j,n,k;

scanf（"%d%f",&n,&a）;

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

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

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

{for（i=k;i<=n;i++）

y[i]=y[k-1]+（y[i]-y[k-1]）*（a-x[k-1]）/（x[i]-x[k-1]）;}

printf（"%f\n",y[n]）;

}

5、复化梯形法

P95.9设,用复化梯形法求积分的近似值

程序：

#include

#include

doublefun（double）;

voidmain（）

{

doublea,b,h,s,x,y;

intn,k;

scanf（"%lf%lf%d",&a,&b,&n）;

h=（b-a）/n;

s=0;

x=a;

y=fun（x）;

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

{s=s+fun（x）;

x=x+h;

s=s+fun（x）;

}

s=（h/2）*s;

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

}

doublefun（doublem）

{doublen;

n=exp（-m）*sin（4*m）+1;

return（n）;

}

运行结果：

6、复化辛甫生算法

P95.9设,用复化辛甫生法求积分的近似值

程序：

#include

#include

doublefun（double）;

voidmain（）

{

doublea,b,h,s,x,y;

intn,k;

scanf（"%lf%lf%d",&a,&b,&n）;

h=（b-a）/n;

s=0;

x=a;

y=fun（x）;

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

{s=s+fun（x）;

x=x+h/2;

s=s+4*fun（x）;

x=x+h/2;

s=s+fun（x）;}

s=（h/6）*s;

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

}

doublefun（doublem）

{doublen;

n=exp（-m）*sin（4*m）+1;

return（n）;

}

运行结果：

7、二阶龙格库塔方法

求解初值问题：

取h=0.2

程序：

#include

#include

floatfun（float,float）;

voidmain（）

{floath,x0,y0,x1,y1,k1,k2;

intn,N;

scanf（"%f%f%f%d",&x0,&y0,&h,&N）;

n=1;

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

{x1=x0+h;

k1=fun（x0,y0）;

k2=fun（x0+h/2,y0+h/2*k1）;

y1=y0+h*k2;

printf（"%f,%f\n",x1,y1）;

x0=x1;y0=y1;}

}

floatfun（floata,floatb）

{floatm;

m=b-2*a/b;

return（m）;

}

运行结果：

8、四阶龙格库塔方法

求解初值问题：

取h=0.2

程序：

#include

#include

floatfun（float,float）;

voidmain（）

{floath,x0,y0,x1,y1,k1,k2,k3,k4;

intn,N;

scanf（"%f%f%f%d",&x0,&y0,&h,&N）;

n=1;

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

{x1=x0+h;

k1=fun（x0,y0）;

k2=fun（x0+h/2,y0+h/2*k1）;

k3=fun（x0+h/2,y0+h/2*k2）;

k4=fun（x1,y0+h*k3）;

y1=y0+h/6*（k1+2*k2+2*k3+k4）;

printf（"%f,%f",x1,y1）;

x0=x1;y0=y1;}

}

floatfun（floata,floatb）

{floatm;

m=b-2*a/b;

return（m）;

}

运行结果：

9、改进的欧拉方法

求解初值问题：

程序：

#include

#include

floatfun（float,float）;

voidmain（）

{floatx0,y0,h,x1,y1,yp,yc;

intn,N;

scanf（"%f%f%f%d",&x0,&y0,&h,&N）;

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

{x1=x0+h;

yp=y0+h*fun（x0,y0）;

yc=y0+h*fun（x1,yp）;

y1=（yp+yc）/2;

printf（"%f,%f",x1,y1）;

x0=x1;

y0=y1;}

}

floatfun（floata,floatb）

{floatm;

m=b-2*a/b;

return（m）;

}

运行结果：

10、迭代法

P131例2用迭代法求方程在附近的一个根，要求精度为

程序：

#include

#include

floatfun（float）;

voidmain（）

{floatx0,x1,c;

 intk,N;

 scanf（"%f%f%d",&x0,&c,&N）;

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

 {x1=fun（x0）;

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

 if（fabs（x1-x0）

 x0=x1;

 }

if（k-1==N）

printf（"Failure!

\n"）;

}

floatfun（floatm）

{floatn;

n=exp（-m）;

return（n）;

}

运行结果：

11、埃特金加速方法：

程序：

#include

#include

floatfun（float）;

voidmain（）

{floatx0,x1,x2,c;

 intk,N;

 scanf（"%f%f%d",&x0,&c,&N）;

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

 {x1=fun（x0）;

 x2=fun（x1）;

 x2=x2-（x2-x1）*（x2-x1）/（x2-2*x1+x0）;

 if（fabs（x2-x0）

 {printf（"%f\n",x2）;

 break;}

 x0=x2;

 }

if（k-1==N）

printf（"Failure!

\n"）;

}

floatfun（floatm）

{floatn;

n=exp（-m）;

return（n）;

}

运行结果：

12、牛顿迭代法：

例5、用牛顿法解方程

牛顿公式为：

，取x=0.5

程序：

#include

#include

floatfun（float）;

floatff（float）;

voidmain（）

{floatx0,x1,c;

 intk,N;

 scanf（"%f%f%d",&x0,&c,&N）;

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

 {if（ff（x0）!

=0）

 {x1=x0-fun（x0）/ff（x0）;

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

 if（fabs（x1-x0）

 x0=x1;}

 else

       printf（"****"）;

 }

if（k==N）

printf（"Failure!

\n"）;

}

floatfun（floatm）

{floatn;

n=m-exp（-m）;

return（n）;

}

floatff（floatm）

{floatn;

n=1+m;

return（n）;

}

运行结果：

13、追赶法：

P198.3.用追赶法求解下列方程组：

程序：

#include

#include

voidmain（）

{floata[100],b[100],c[100],d[100],t;

inti,n;

scanf（"%d",&n）;

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

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

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

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

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

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

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

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

c[1]=c[1]/b[1];

d[1]=d[1]/b[1];

for（i=2;i

{t=b[i]-c[i-1]*a[i];

c[i]=c[i]/t;

d[i]=（d[i]-d[i-1]*a[i]）/t;}

d[n]=（d[n]-d[n-1]*a[n]）/（b[n]-c[n-1]*a[n]）;

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

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

{d[i]=d[i]-c[i]*d[i+1];

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

}

运行结果：

14、雅克比迭代

程序：

#include

#include

#defineN50                             

#defineM4

voidmain（）

{doublex[M],y[M],a[M][M],b[M],d[M],c,t;

 doubleff（double[],int）;

 intk,i,j,n;

n=M-1;

scanf（"%lf",&c）;                                     

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

{scanf（"%lf",&x[i]）;

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

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

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

scanf（"%lf",&a[0][i*M+j]）;

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

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

{for（j=1,t=0;j<=n;j++）

{if（j==i）

continue;

else

t=t+a[i][j]*x[j];}

y[i]=（b[i]-t）/a[i][i];

}

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

d[i]=fabs（x[i]-y[i]）;

t=ff（d,n+1）;

if（t

else

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

x[i]=y[i];

}

if（k==N）

printf（"Failure!

\n"）;

if（k

{printf（"k=%d\n",k）;

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

printf（"y[i]=%f\n",y[i]）;}

}

doubleff（doublea[],intn）

 {doublep;

 intt;

 p=a[1];

 for（t=2;t

{if（p

p=a[t];}

return（p）;

}

运行结果：

15、蛋白质设计：

程序：

#include

#include

#include

 voidmain（）

{printf（"横坐标    纵坐标   竖坐标\n"）;  

FILE*fp;

charc[100],x[3000][7],y[3000][7],z[3000][7];

floata[3000],b[3000],d[3000],r[3000];

inti,k=0,j,n=0,m,p;

floatX[3000],Y[3000],Z[3000],s1=0,s2=0,s3=0,av_x,av_y,av_z;

fp=fopen（"1a1c.pdb","r"）;

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

{

fgets（c,81,fp）;

if（c[0]=='A'&&c[1]=='T'&&c[2]=='O'&&c[3]=='M'）  

{

for（j=0;j<6;j++）

{x[k][j]=c[32+j];                               

printf（"%c",x[k][j]）;

}

printf（"   "）;

for（j=0;j<6;j++）                                

{y[k][j]=c[40+j];                               

printf（"%c",y[k][j]）;}

printf（"   "）;

for（j=0;j<6;j++）                               

{z[k][j]=c[48+j];

printf（"%c",z[k][j]）;

}

printf（"   "）;

X[k]=atof（x[k]）;

s1+=X[k];

Y[k]=atof（y[k]）;

s2+=Y[k];

Z[k]=atof（z[k]）;

s3+=Z[k];

printf（"\n"）;

k++;

if（c[77]=='C'）n++;

}

}

av_x=s1/k;

av_y=s2/k;

av_z=s3/k;

printf（"共有:

%d\n",k+1）;

printf（"av_x=%f,av_y=%f,av_z=%f\n",av_x,av_y,av_z）;

printf（"C原子个数为:

%d\n",n）;

printf（"平移原点后的坐标:

\n"）;

for（m=0;m

{a[m]=X[m]-av_x;printf（"%f    ",a[m]）;              

 b[m]=Y[m]-av_y;printf（"%f    ",b[m]）;

 d[m]=Z[m]-av_z;printf（"%f    ",d[m]）;

 r[m]=sqrt（a[m]*a[m]+b[m]*b[m]+d[m]*d[m]）;printf（"%f    ",r[m]）;

printf（"\n"）;

}

}

17高斯消去法:

#include"stdio.h"

#include"math.h"

main（）

{doublea[3][3]={1,1,1,0,4,-1,2,-2,1},b[3]={6,5,1},x[10]={0};

 inti,j,k,n=3;

 for（k=0;k

 { for（i=k+1;i

   { for（j=k+1;j

     { a[i][j]=a[i][j]-a[k][j]*a[i][k]/a[k][k];

     }

     b[i]=b[i]-b[k]*a[i][k]/a[k][k];

   }

 }

 x[n-1]=b[n-1]/a[n-1][n-1];

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

 { k=n-i; 

 for（j=k+1;j

   { x[k]+=a[k][j]*x[j];

   }

   x[k]=（b[k]-x[k]）/a[k][k];

 }

for（k=0;k

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

}