一、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 ;}