文档视界
文档视界 最新最全的文档下载
当前位置：文档视界 › 数值计算方法实验程序

数值计算方法实验程序

一、Newton插值法

#include

#define MAX_N 20

typedef struct tagPOINT

{double x;

double y;

}POINT;



int main()

{int n,i,j;

POINT points[MAX_N+1];double diff[MAX_N+1];

double x,tmp,newton=0;

printf("\nInput n value:");

scanf("%d",&n);

printf("Now input the (x_i,y_i),....,%d:\n",n);

for(i=0;i<=n;i++)

scanf("%lf%lf",&points[i].x,&points[i].y);

printf("Now input the x value:");

scanf("%lf",&x);

for(i=0;i<=n;i++) diff[i]=points[i].y;

for(i=0;i<=n;i++)

{for(j=n;j>i;j--)

{diff[j]=(diff[j]-diff[j-1])/(points[j].x-points[j-1-i].x);

}

}

tmp=1;newton=diff[0];

for(i=0;i
{tmp=tmp*(x-points[i].x);

newton=newton+tmp*diff[i+1];

}

printf("newton(%f)=%f\n",x,newton);

return 0;

}

二、Lagrange插值法

double Lagrange(int n,double u,double x[],double y[])

{int i,j;

double l,Ln=0;

for(i=0;i<=n;i++)

{l=1.0;

for(j=0;j<=n;j++)

if (i!=j)l=l*(u-x[j])/(x(i)-x[j]);

Ln=Ln+l*y[i];

}

return(Ln);

}

#include

#define N 20

Void main(void)

{int i,n;

double x[N],y[N],u.fu;

double Lagrange(int n,double u,double x[],double y[]);

printf(“请输入插值多项式的次数N及插入点U：\n”);

scanf(“%d,%lf”,&n,&u);

printf(“请输入n+1个插值节点x,y:\n”);

for(i=0;i<=n;i++)

scanf(“%lf,%lf”,&x[i],&y[i]);

fu=Lagrange(n,u,x,y);

printf(“插入点u的近似值为%lf\n”,fu);



}

三、用梯形公式

#include

#include

double f(double x)

{return(sqrt(pow(x,4)+1));

}

double T(double a,double b)

{double u,h;

h=(b-a)/2;

u=f(a)+f(b);

return (h*u);

}

void main()

{ double a,b;

printf("请输入积分区域：a,b\n");

scanf("%lf,%lf",&a,&b);

double T(double ,double );

printf("%lf,%lf",T(a,b));

四、Simpson公式求积分

#include

#include

double f(double x)

{

return(sqrt(2-pow(sin(x),2)));

}

double S(double a,double b)

{

double u,h,c;

h=(b-a)/6;

c=(a+b)/2;

u=f(a)+f(b)+4*f(c);

return (h*u);

}

void main()

{ double a,b;

printf("请输入积分区域：a,b\n");

scanf("%lf,%lf",&a,&b);

double S(double ,double );

printf("%lf ",S(a,b))







五、复化数值积分公式的自动控制误差算法



复化Simpson自动控制误差

#include

double AS(double (*f)(double x),double a,double b,double e)

{int i;int n=2;

double h,Jn,J2n,Sn,S2n;

h=(b-a)/n;

Jn=(*f)((a+b)/2);

Sn=h/3*((*f)(a)+(*f)(b)+4*Jn);

while(1)

{J2n=0;

for(i=1;i<=n;i++)

J2n=J2n+(*f)(a+(2*i-1)*h/2);

S2n=Sn/2+h/6*(4*J2n-2*Jn);

if(fabs(S2n-Sn)<15*e){return S2n;break;}

else{ h=h/2;n=2*n;J2n=Jn;Sn=S2n;}

}}

复化梯形公式

double comTx(double (*f)(double),double a,double b,int n)

{int i; double h,Tn;

h=(b-a)/n;

Tn=(*f)(a)+(*f)(b);


for(i=1;i<=n-1;i++)

Tn=Tn+2*(*f)(a+i*h);

Tn=(h/2)*Tn;

return(Tn);}

复化梯形自动控制误差

double AT(double (*f)(double x),double a,double b,double e)

{int i; int n=10;

double h,Js,Tn,T2n;

h=(b-a)/n;

Tn=comTx(f,a,b,n);

while(1)

{Js=0;

for(i=1;i<=n;i++)

Js=Js+(*f)(a+(2*i-1)*h/2);

T2n=Tn/2+h/2*Js;

if(fabs(T2n-Tn)<3*e)break;

else{ h=h/2;n=2*n;Tn=T2n;}

}

return T2n;

}

六、Romberg公式计算积分

#include "stdio.h"

#include "math.h"

#define MAX 20

double f(double x)//要求的积分公式

{ return(lnx);

}

double computeTn(double (*f)(double x),double a,double b, int n)//复化梯形公式

{int i;

double sum=0;

double h=(b-a)/n;

for(i=1;i
sum+=f(a+i*h);

sum+=(f(a)+f(b))/2;

return (h*sum);

}

double Romberg(double (*f)(double x),double a,double b,double ep)//Romberg求积分

{ int j,n=1,k;

double R[3][MAX];

R[2][1]=computeTn(f,a,b,n);

printf("%lf\n",R[2][1]);

for(k=2;k
{ for(j=1;j
R[1][j]=R[2][j];

n*=2;

R[2][1]=computeTn(f,a,b,n);

printf("%lf\t",R[2][1]);

for(j=2;j<=k;j++)

{ R[2][j]=R[2][j-1]+(R[2][j-1]-R[1][j-1])/(pow(4,j-1)-1);

printf("%lf\t",R[2][j]);

}

printf("\n");

if(fabs(R[1][k-1]-R[2][k])
}

k--;

if(fabs(R[1][k-1]-R[2][k])>ep)

{printf("Return no solved...\n");

return 0;

}

}

七、用牛顿迭代法、弦截法求非线性方程的根

///用牛顿迭代法的C源程序///

#include

#include

#define m 20

double f(double x)

{return((x*(x*(x-7.7)+19.2)-15.3));}

double f1(double x)

{return(x*(3*x-15.4)+19.2);}

main()

{int k;

double x,x0,e;

printf("请输入x0和允许误差e的值:\n");

scanf("%lf,%lf",&x0,&e);

for(k=1;k<=m;k++)

{x=x0-f(x0)/f1(x0);

if(fabs(x-x0)
x0=x;

}

printf("给定的迭代次数太小.停机");

}

///用弦截法的C源程序///

#include

#include

#define m 20

double f(double x)

{return((x*(x*(x-7.7)+19.2)-15.3));}

main()

{int k;

double x,x0,x1,f0,f1,e;

printf("请输入x0,x1和允许误差e的值:\n");

scanf("%lf,%lf,%lf",&x0,&x1,&e);

f0=f(x0);

for(k=1;k<=m;k++)

{f1=f(x1);

x=x1-(f1*(x1-x0)/(f1-f0));

if(fabs(x-x1)
f0=f1;

x0=x1;

x1=x;}

printf("给定的迭代次数太小.停机");

}

八、Gauss-Seidel迭代法解线性方程组的算法

#include

#include

#define eps 0.000001

#define max 100

#define N 20

double norm_inf(double x[],int n)

{ double norm;

int i;

norm=fabs(x[0]);

for(i=1;i
{

if(fabs(x[i])>norm)

norm=fabs(x[i]);

}

return norm;

}

void seidel(double a[N][N],double g[N],int n)

{

double b[N][N]={0},x0[

N],x1[N],x1_x0[N],norm,temp;

int i,j,k;

for(i=0;i
{

g[i]=g[i]/a[i][i];

for(j=0;j
{ if(i==j)

continue;

b[i][j]=-a[i][j]/a[i][i];

}

}

for(i=0;i
{

x0[i]=0;

x1[i]=1;

x1_x0[i]=x1[i]-x0[i];

}

k=0;

norm=norm_inf(x1_x0,n);

while((norm>=eps)&&(k
{

for(i=0;i
x0[i]=x1[i];

for(i=0;i
{

temp=0;

for(j=0;j<=i-1;j++)

temp=temp+b[i][j]*x1[j];

for(j=i+1;j
temp=temp+b[i][j]*x0[j];

x1[i]=temp+g[i];

x1_x0[i]=x1[i]-x0[i];

}

norm=norm_inf(x1_x0,n);

k++;

}

for(i=0;i
printf("x[%d]=%lf\n",i,x1[i]);

printf("%d times iteration.\n",k);

}

int main()

{

double b[N][N]={{10,-2,-1},{-2,10,-1},{-1,-2,5}},f[N]={0,-21,-20};

printf("\n");

seidel(b,f,3);

return 0;

}

九、用列主元消元法

#include

#include

#define n 3

void print(double a[n][n+1]);

void gauss(double a[n][n+1],double x[n])

{ int i,j,k;

double temp,s,l;

for(i=0;i
{ k=i;

for(j=i+1;j
{ if(fabs(a[j][i])>fabs(a[k][i]))

k=j;

}

if(k!=i)

for(j=i;j<=n;j++)

{temp=a[i][j];

a[i][j]=a[k][j];

a[k][j]=temp;

}

for(j=i+1;j
{ l=1.0*a[j][i]/a[i][i];

for(k=0;k
a[j][k]=a[j][k]-a[i][k]*l;

}

print(a);

printf("lf");}

print(a);

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

for(i=n-2;i>=0;i--)

{ s=0.0;

for(j=i;j
{ if(j==i)

continue;

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

x[i]=(a[i][n]-s)/a[i][i];

}}

void print(double a[n][n+1])

{int i,j;

for(i=0;i
{ for(j=0;j
printf("%lf",a[i][j]);

printf("lf");}}

十、Doolittle分解法

#include "iostream.h"

#include "stdio.h"

#define max 20

void main()

{int i,k,n,j,p;

double x,v,y;

double a[max][max],r[max][max],l[max][max];

double b[max],Y[max],X[max];

printf("Input n:"); scanf("%d",&n);

printf("Input a[%d][%d]:\n",n,n);

for(i=1;i
for(j=1;j
{ scanf("%lf",&x); a[i][j]=x; }

printf("b[i]:");

for(i=1;i
{ scanf("%lf",&y); b[i]=y; }

for(k=1;k
{ for(i=k;i
{ v=0;

for(p=1;p
v=l[k][p]*r[p][i]+v;

r[k][i]=a[k][i]-v; printf("r[%d][%d]=%f\n",k,i,r[k][i]); }

for(i=k+1;i
{ v=0;

for(p=1;p
v=v+l[i][p]*r[p][k];

l[i][k]=(a[i][k]-v)/r[k][k]; printf("l[%d][%d]=%f\n",i,k,l[i][k]);}

for(i=1;i
{ v=0;

for(k=1;k
v=v+l[i][k]*Y[k];

Y[i]=b[i]-v;}

for(i=n;i>0;i--)

{ v=0;

for(k=i+1;k
v=v+r[i][k]*X[k];

X[i]=(Y[i]-v)/r[i][i];}

printf("\n");

for(i=1;i


printf("x[%d]=%f\n",i,X[i]);

printf("\n");

return ;}


相关文档
相关文档 最新文档

侵权投诉  © 2013-2023 www.docsj.com  站点地图 | 手机版

闽ICP备11023808号-7  本站文档均来自互联网及网友上传分享，本站只负责收集和整理，有任何问题可通过上访投诉通道进行反馈