数值分析源码和库

编前语：

这次放出的是只有源码，大部分是拼凑来的，我都做了修改（郁闷那么多错还发到网上共享），多为通用算法，专业级的算法有，乘我还是学生我会想办法为大家弄的。以后会放出lib的供大家连接，当然lib版的更稳定些，更全面些，还会加入各种异常，便于大家调试。好了，不说了。

// JKs_NumericalMethod.h:
//
//////////////////////////////////////////////////////////////////////

#if !defined(__JKs_NumericalMethod_H_)
#define __JKs_NumericalMethod_H_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000

namespace jks
{
//////////////////////////////////////////////////////////////////////////
/// 二分法
extern double EPS_Bisection; //根容许的误差
extern double DELTA_Bisection; //|f(x)|容许的误差
float Bisection( float a,float b,float (*f)(float) );
double Bisection( double a,double b,double (*f)(double) );

///复化辛卜生公式
float Simpson(float (*f)(float),float a,float b,int n);
double Simpson(double (*f)(double),double a,double b,int n);

//改进欧拉法
void ModEuler(float (*f1)(float,float),float x0,float y0,float xn,int n);
void ModEuler(double (*f1)(double,double),double x0,double y0,double xn,int n);

//高斯-赛德尔迭代法
extern double N_GauseSeidel;
float *GauseSeidel(float *a,int n);
double *GauseSeidel(double *a,int n);

//拉格郎日插值多项式
float Lagrange(float *x,float *y,float xx,int n);
double Lagrange(double *x,double *y,double xx,int n);

//列主元高斯消去法
float *ColPivot( float *c,int n );
double *ColPivot( double *c,int n );

//龙贝格算法
float Romberg(float a,float b,float (*f)(float),float epsilon);
double Romberg(double a,double b,double (*f)(double),double epsilon);

//龙格-库塔算法
void Runge_Kutta(float (*f)(float x,float y),float a,float b,float y0,int N);
void Runge_Kutta(double (*f)(double x,double y),double a,double b,double y0,int N);

//幂法
//extern double N_PowerMethod;
extern double EPS_PowerMethod;
extern double KM_PowerMethod;
void PowerMethod(float *A);
void PowerMethod(double *A);

//牛顿迭代法
extern double N_Newton;
extern double EPS_Newton;
extern double ETA_Newton;
float Newton(float (*f)(float),float (*f1)(float),float x0);
double Newton(double (*f)(double),double (*f1)(double),double x0);

//牛顿值多项式
//extern double N_Difference;
void Difference(float *x,float *y,int n);
void Difference(double *x,double *y,int n);

//四阶阿当姆斯预测-校正公式
void Adams(float a,float b,int N,float (*f)(float x,float y),float y0);
void Adams(double a,double b,int N,double (*f)(double x,double y),double y0);

//雅可比迭代法
extern double EPS_Jacobi;
extern double MAX_Jacobi ;
float *Jacobi(float a[3][4],int n);
double *Jacobi(double a[3][4],int n);

//自适应梯形公式（变步长）
float AutoTrap(float (*f)(float),float a,float b);
double AutoTrap(double (*f)(double),double a,double b);

//最小二乘法
float *Approx(float *x,float *y,int m,int n);
double *Approx(double *x,double *y,int m,int n);
//////////////////////////////////////////////////////////////////////////
} //
#endif // !defined(__JKs_NumericalMethod_H_

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

// JKs_NumericalMethod.cpp:
//
//////////////////////////////////////////////////////////////////////

#include<iostream>
#include<math.h>
using namespace std;
#include "JKs_NumericalMethod.h"

/*
二分法
复化辛卜生公式
改进欧拉法
高斯-赛德尔迭代法
拉格郎日插值多项式
列主元高斯消去法
龙贝格算法
龙格-库塔算法
幂法
牛顿迭代法
牛顿值多项式
四阶阿当姆斯预测-校正公式
雅可比迭代法
自适应梯形公式（变步长）
最小二乘法
*/
namespace jks
{
//////////////////////////////////////////////////////////////////////////
/// 二分法
double EPS_Bisection =5e-6; //根容许的误差
double DELTA_Bisection =1e-6; //|f(x)|容许的误差
float Bisection( float a,float b,float (*f)(float) )
{
float c,fc,fa=f(a),fb=f(b);
int n=1;
cout<<"二分次数/tc/tf(c)/n";
while(1)
{
if( fa*fb>0 )
return 0;
c =(a+b)/2,fc =f(c);
if( fabs(fc) < DELTA_Bisection )
break;
else if (fa*fc<0)
{
b=c;
fb=fc;
}
else
{
a=c;
fa=fc;
}
if( b-a< EPS_Bisection )
break;
cout<<'/t'<<n++<<'/t'<<c<<'/t'<<fc<<endl;
}
return c;
}

float f_Bisection(float x)
{
return x*x*x+x*x-3*x-3; //求解的方程
}

//
double Bisection( double a,double b,double (*f)(double) )
{
double c,fc,fa=f(a),fb=f(b);
int n=1;
cout<<"二分次数/tc/tf(c)/n";
while(1)
{
if( fa*fb>0 )
return 0;
c =(a+b)/2,fc =f(c);
if( fabs(fc) < DELTA_Bisection )
break;
else if (fa*fc<0)
{
b=c;
fb=fc;
}
else
{
a=c;
fa=fc;
}
if( b-a< EPS_Bisection )
break;
cout<<'/t'<<n++<<'/t'<<c<<'/t'<<fc<<endl;
}
return c;
}

double f_Bisection(double x)
{
return x*x*x+x*x-3*x-3; //求解的方程
}

//////////////////////////////////////////////////////////////////////////
///复化辛卜生公式
float Simpson(float (*f)(float),float a,float b,int n)
{
int k;
float s,s1,s2=0.0;
float h=(b-a)/n;
s1=f(a+h/2);
for(k=1;k<=n-1;k++)
{
s1+=f(a+k*h+h/2);
s2+=f(a+k*h);
}
s=h/6*(f(a)+4*s1+2*s2+f(b));
return s;
}
float f_Simpson(float x)
{
return 1/(1+x*x);
/* if(x==0)
return 1;
else
return sin(x)/x;*/
}

//
double Simpson(double (*f)(double),double a,double b,int n)
{
int k;
double s,s1,s2=0.0;
double h=(b-a)/n;
s1=f(a+h/2);
for(k=1;k<=n-1;k++)
{
s1+=f(a+k*h+h/2);
s2+=f(a+k*h);
}
s=h/6*(f(a)+4*s1+2*s2+f(b));
return s;
}
double f_Simpson(double x)
{
return 1/(1+x*x);
/* if(x==0)
return 1;
else
return sin(x)/x;*/
}

//////////////////
void test_Simpson()
{
int i,n=2;
float s;
float f_Simpson(float);
float Simpson(float (*)(float),float,float,int );
for( i=0;i<=2;i++ )
{
s=Simpson(f_Simpson,0,1,n);
cout<<"s("<<n<<")="<<s<<endl;
n*=2;
}
}

void test_Simpson2()
{
int i,n=2;
double s;
double f_Simpson(double);
double Simpson(double (*)(double),double,double,int );
for( i=0;i<=2;i++ )
{
s=Simpson(f_Simpson,0,1,n);
cout<<"s("<<n<<")="<<s<<endl;
n*=2;
}
}
//////////////////////////////////////////////////////////////////////////
//改进欧拉法
void ModEuler(float (*f1)(float,float),float x0,float y0,float xn,int n)
{
int i;
float yp,yc,x=x0,y=y0,h=(xn-x0)/n;
cout<<"x[0]="<<x<<'/t'<<"y[0]"<<y<<endl;
for(i=1;i<=n;i++)
{
yp=y+h*f1(x,y);
x=x0+i*h;
yc=y+h*f1(x,yp);
y=(yp+yc)/2.0;
cout<<"x["<<i<<"]="<<x<<" y["<<i<<"]="<<y<<endl;
}
}

void ModEuler(double (*f1)(double,double),double x0,double y0,double xn,int n)
{
int i;
double yp,yc,x=x0,y=y0,h=(xn-x0)/n;
cout<<"x[0]="<<x<<'/t'<<"y[0]"<<y<<endl;
for(i=1;i<=n;i++)
{
yp=y+h*f1(x,y);
x=x0+i*h;
yc=y+h*f1(x,yp);
y=(yp+yc)/2.0;
cout<<"x["<<i<<"]="<<x<<" y["<<i<<"]="<<y<<endl;
}
}

double N_ModEuler =20;
float f1_ModEuler(float x,float y)
{
return -x*y*y;
}
void test_ModEuler()
{

float xn=5.0,x0=0.0,y0=2.0;
float f1_ModEuler(float ,float);
ModEuler(f1_ModEuler,x0,y0,xn,N_ModEuler);
}

double f1_ModEuler(double x,double y)
{
return -x*y*y;
}
void test_ModEuler2()
{

double xn=5.0,x0=0.0,y0=2.0;
double f1_ModEuler(double ,double);
ModEuler(f1_ModEuler,x0,y0,xn,N_ModEuler);
}

//////////////////////////////////////////////////////////////////////////
//高斯-赛德尔迭代法
double N_GauseSeidel =100;
float *GauseSeidel(float *a,int n)
{

int i,j,nu=0;
float *x,dx;
x =new float
;
for( i=0;i<=n-1;i++ )
x[i]=0.0;
do
{
for( i=0;i<=n-1;i++ )
{
dx=0.0;
for( j=0;j<=n-1;j++ )
dx+=*(a+i*(n+1)+j)*x[j];
dx=(*(a+i*(n+1)+n)-dx)/(*(a+i*(n+1)+i));
x[i]+=dx;
}
if( nu>N_GauseSeidel )
{
cout<<"迭代发散/n";
return 0;
}
nu++;
}
while(fabs(dx)>1e-6);
return x;
}

double *GauseSeidel(double *a,int n)
{

int i,j,nu=0;
double *x,dx;
x =new double
;
for( i=0;i<=n-1;i++ )
x[i]=0.0;
do
{
for( i=0;i<=n-1;i++ )
{
dx=0.0;
for( j=0;j<=n-1;j++ )
dx+=*(a+i*(n+1)+j)*x[j];
dx=(*(a+i*(n+1)+n)-dx)/(*(a+i*(n+1)+i));
x[i]+=dx;
}
if( nu>N_GauseSeidel )
{
cout<<"迭代发散/n";
return 0;
}
nu++;
}
while(fabs(dx)>1e-6);
return x;
}

void test_GauseSeidel()
{
int i;
float *x;
float c[12]={8.0,-3.0,2.0,20.0,
4.0,11.0,-1.0,33.0,
6.0,3.0,12.0,36.0}; //here must be modifed
float *GauseSeidel(float *,int);
x=GauseSeidel(c,3);
for( i=0;i<=2;i++ )
cout<<"x["<<i<<"]="<<x[i]<<endl;
}

void test_GauseSeidel2()
{
int i;
double *x;
double c[12]={8.0,-3.0,2.0,20.0,
4.0,11.0,-1.0,33.0,
6.0,3.0,12.0,36.0}; //here must be modifed
double *GauseSeidel(double *,int);
x=GauseSeidel(c,3);
for( i=0;i<=2;i++ )
cout<<"x["<<i<<"]="<<x[i]<<endl;
}

//////////////////////////////////////////////////////////////////////////
//拉格郎日插值多项式
float Lagrange(float *x,float *y,float xx,int n)
{
int i,j;
float *a,yy=0.0;
a =new float
;
for( i=0;i<=n-1;i++ )
{
a[i]=y[i];
for( j=0;j<=n-1;j++ )
if(j!=i)
a[i]*=(xx-x[j])/(x[i]-x[j]);
yy+=a[i];
}
delete[] a;
return yy;
}

double Lagrange(double *x,double *y,double xx,int n)
{
int i,j;
double *a,yy=0.0;
a =new double
;
for( i=0;i<=n-1;i++ )
{
a[i]=y[i];
for( j=0;j<=n-1;j++ )
if(j!=i)
a[i]*=(xx-x[j])/(x[i]-x[j]);
yy+=a[i];
}
delete[] a;
return yy;
}

void test_Lagrange()
{
float x[4]={0.56160,0.56280,0.56401,0.56521};
float y[4]={0.82741,0.82659,0.82577,0.82495};
float xx=0.5635,yy;
float Lagrange(float *,float *,float ,int);
yy =Lagrange(x,y,xx,4);
cout<<"x="<<xx<<'/t'<<"y="<<yy<<endl;
}

void test_Lagrange2()
{
double x[4]={0.56160,0.56280,0.56401,0.56521};
double y[4]={0.82741,0.82659,0.82577,0.82495};
double xx=0.5635,yy;
double Lagrange(double *,double *,double ,int);
yy =Lagrange(x,y,xx,4);
cout<<"x="<<xx<<'/t'<<"y="<<yy<<endl;
}

//////////////////////////////////////////////////////////////////////////
//列主元高斯消去法
float *ColPivot( float *c,int n )
{
int i,j,t,k;
float *x,p;
x=new float
;
for( i=0;i<=n-2;i++)
{
k=i;
for(j=i+1;j<=n-1;j++)
if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))
k=j;
if(k!=i)
for( j=i;j<=n;j++ )
{
p=*(c+i*(n+1)+j);
*(c+i*(n+1)+j)=*(c+k*(n+1)+j);
*(c+k*(n+1)+j)=p;
}
for( j=i+1;j<=n-1;j++ )
{
p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));
for( t=i;t<=n;t++ )
*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));
}
}
for( i=n-1;i>=0;i--)
{
for( j=n-1;j>=i+1;j--)
(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));
x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));
}
return x;
}

double *ColPivot( double *c,int n )
{
int i,j,t,k;
double *x,p;
x=new double
;
for( i=0;i<=n-2;i++)
{
k=i;
for(j=i+1;j<=n-1;j++)
if(fabs(*(c+j*(n+1)+i))>(fabs(*(c+k*(n+1)+i))))
k=j;
if(k!=i)
for( j=i;j<=n;j++ )
{
p=*(c+i*(n+1)+j);
*(c+i*(n+1)+j)=*(c+k*(n+1)+j);
*(c+k*(n+1)+j)=p;
}
for( j=i+1;j<=n-1;j++ )
{
p=(*(c+j*(n+1)+i))/(*(c+i*(n+1)+i));
for( t=i;t<=n;t++ )
*(c+j*(n+1)+t)-=p*(*(c+i*(n+1)+t));
}
}
for( i=n-1;i>=0;i--)
{
for( j=n-1;j>=i+1;j--)
(*(c+i*(n+1)+n))-=x[j]*(*(c+i*(n+1)+j));
x[i]=*(c+i*(n+1)+n)/(*(c+i*(n+1)+i));
}
return x;
}



void test_ColPivot()
{
int i;
float *x;
float c[3][4] ={0.101,2.304,3.555,1.183,
-1.347,3.712,4.623,2.137,
-2.835,1.072,5.643,3.035};
float *ColPivot(float *,int);
x=ColPivot(c[0],3);
for( i=0;i<=2;i++ )
cout<<"x("<<i<<")="<<x[i]<<endl;
}

void test_ColPivot2()
{
int i;
double *x;
double c[3][4] ={0.101,2.304,3.555,1.183,
-1.347,3.712,4.623,2.137,
-2.835,1.072,5.643,3.035};
double *ColPivot(double *,int);
x=ColPivot(c[0],3);
for( i=0;i<=2;i++ )
cout<<"x("<<i<<")="<<x[i]<<endl;
}

//////////////////////////////////////////////////////////////////////////
//龙贝格算法
float Romberg(float a,float b,float (*f)(float),float epsilon)
{
int n=1,k;
float h=b-a,x,temp;
float T1,T2,S1,S2,C1,C2,R1,R2;
T1=(b-a)/2*((*f)(a)+(*f)(b));
while(1)
{
temp=0;
for(k=0;k<=n-1;k++)
{
x=a+k*h+h/2;
temp+=(*f)(x);
}
T2=(T1+temp*h)/2;
if(fabs(T2-T1)<epsilon)
return T2;
S2=T2+(T2-T1)/3.0;
if(n==1)
{
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
C2=S2+(S2-S1)/15;
if(n==2)
{
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
R2=C2+(C2-C1)/63;
if(n==4)
{
R1=R2;
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
if(fabs(R2-R1)<epsilon)
return R2;
R1=R2;
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
}
}

double Romberg(double a,double b,double (*f)(double),double epsilon)
{
int n=1,k;
double h=b-a,x,temp;
double T1,T2,S1,S2,C1,C2,R1,R2;
T1=(b-a)/2*((*f)(a)+(*f)(b));
while(1)
{
temp=0;
for(k=0;k<=n-1;k++)
{
x=a+k*h+h/2;
temp+=(*f)(x);
}
T2=(T1+temp*h)/2;
if(fabs(T2-T1)<epsilon)
return T2;
S2=T2+(T2-T1)/3.0;
if(n==1)
{
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
C2=S2+(S2-S1)/15;
if(n==2)
{
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
R2=C2+(C2-C1)/63;
if(n==4)
{
R1=R2;
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
continue;
}
if(fabs(R2-R1)<epsilon)
return R2;
R1=R2;
C1=C2;
T1=T2;
S1=S2;
h/=2;
n*=2;
}
}

float f_Romberg(float x)
{
return 1/(1+x*x);
}

double f_Romberg(double x)
{
return 1/(1+x*x);
}

void test_Romberg()
{
float epsilon=5e-6;
cout<<"R="<<Romberg(0,1,f_Romberg,epsilon)<<endl;
}

void test_Romberg2()
{
double epsilon=5e-6;
cout<<"R="<<Romberg(0,1,f_Romberg,epsilon)<<endl;
}

//////////////////////////////////////////////////////////////////////////
//龙格-库塔算法
#define N_PowerMethod 3
double EPS_PowerMethod =1e-6;
double KM_PowerMethod =30;
void Runge_Kutta(float (*f)(float x,float y),float a,float b,float y0,int N)
{
float x=a,y=y0,K1,K2,K3,K4;
float h=(b-a)/N;
int i;
cout<<"x[0]="<<x<<'/t'<<"y[0]="<<y<<endl;
for(i=1;i<=N;i++)
{
K1=f(x,y);
K2=f(x+h/2,y+h*K1/2);
K3=f(x+h/2,y+h*K2/2);
K4=f(x+h,y+h*K3);
y=y+h*(K1+2*K2+2*K3+K4)/6;
x=a+i*h;
cout<<"x["<<i<<"]="<<x<<" y["<<i<<"]="<<y<<endl;
}
}

void Runge_Kutta(double (*f)(double x,double y),double a,double b,double y0,int N)
{
double x=a,y=y0,K1,K2,K3,K4;
double h=(b-a)/N;
int i;
cout<<"x[0]="<<x<<'/t'<<"y[0]="<<y<<endl;
for(i=1;i<=N;i++)
{
K1=f(x,y);
K2=f(x+h/2,y+h*K1/2);
K3=f(x+h/2,y+h*K2/2);
K4=f(x+h,y+h*K3);
y=y+h*(K1+2*K2+2*K3+K4)/6;
x=a+i*h;
cout<<"x["<<i<<"]="<<x<<" y["<<i<<"]="<<y<<endl;
}
}

float f_Runge_Kutta(float x,float y)
{
return -x*y*y;
}

double f_Runge_Kutta(double x,double y)
{
return -x*y*y;
}

void test_Runge_Kutta()
{
void Runge_Kutta(float (*f)(float,float),float a,float b,float y0,int N);
float f_Runge_Kutta(float,float);
float a=0,b=5,y0=2;
Runge_Kutta(f_Runge_Kutta,a,b,y0,20);
}

void test_Runge_Kutta2()
{
void Runge_Kutta(double (*f)(double,double),double a,double b,double y0,int N);
double f_Runge_Kutta(double,double);
double a=0,b=5,y0=2;
Runge_Kutta(f_Runge_Kutta,a,b,y0,20);
}

//////////////////////////////////////////////////////////////////////////
//幂法
float MaxValue(float *x,int n)
{
float Max=x[0];
int i;
for (i=1;i<n;i++)
if(fabs(x[i])>fabs(Max))
Max=x[i];
return Max;
}
double MaxValue(double *x,int n)
{
double Max=x[0];
int i;
for (i=1;i<n;i++)
if(fabs(x[i])>fabs(Max))
Max=x[i];
return Max;
}
void PowerMethod(float *A)
{
float MaxValue(float *,int );
float U[N_PowerMethod],V[N_PowerMethod],r2,r1;
float temp;
int i,j,k=0;
for (i=0;i<N_PowerMethod;i++)
U[i]=1;
while(k<KM_PowerMethod)
{
k++;
for(i=0;i<N_PowerMethod;i++)
{
temp=0;
for(j=0;j<N_PowerMethod;j++)
temp+=*(A+i*N_PowerMethod+j)*U[j];
V[i]=temp;
}
for (i=0;i<N_PowerMethod;i++)
U[i]=V[i]/MaxValue(V,N_PowerMethod);
if(k==1)
{ //书上程序有bug，这里就加入下面两句
r1=MaxValue(V,N_PowerMethod);
continue;
}
r2=MaxValue(V,N_PowerMethod);
if(fabs(r2-r1)<EPS_PowerMethod)
break;
r1=r2;
}
cout<<"/nr="<<r2<<endl;
for (i=0;i<N_PowerMethod;i++)
cout<<"/nx["<<i+1<<"]="<<U[i];
}
void PowerMethod(double *A)
{
double MaxValue(double *,int );
double U[N_PowerMethod],V[N_PowerMethod],r2,r1;
double temp;
int i,j,k=0;
for (i=0;i<N_PowerMethod;i++)
U[i]=1;
while(k<KM_PowerMethod)
{
k++;
for(i=0;i<N_PowerMethod;i++)
{
temp=0;
for(j=0;j<N_PowerMethod;j++)
temp+=*(A+i*N_PowerMethod+j)*U[j];
V[i]=temp;
}
for (i=0;i<N_PowerMethod;i++)
U[i]=V[i]/MaxValue(V,N_PowerMethod);
if(k==1)
{ //书上程序有bug，这里就加入下面两句
r1=MaxValue(V,N_PowerMethod);
continue;
}
r2=MaxValue(V,N_PowerMethod);
if(fabs(r2-r1)<EPS_PowerMethod)
break;
r1=r2;
}
cout<<"/nr="<<r2<<endl;
for (i=0;i<N_PowerMethod;i++)
cout<<"/nx["<<i+1<<"]="<<U[i];
}

void test_PowerMethod()
{
float A[N_PowerMethod][N_PowerMethod]={{2,-1,0},{-1,2,-1},{0,-1,2}};
PowerMethod(A[0]);
}

void test_PowerMethod2()
{
double A[N_PowerMethod][N_PowerMethod]={{2,-1,0},{-1,2,-1},{0,-1,2}};
PowerMethod(A[0]);
}

//////////////////////////////////////////////////////////////////////////
//牛顿迭代法
double N_Newton =100;
double EPS_Newton =1e-6;
double ETA_Newton =1e-8;
float Newton(float (*f)(float),float (*f1)(float),float x0)
{
float x1,d;
int k=0;
do{
x1=x0-f(x0)/f1(x0);
if((k++>N_Newton)||(fabs(f1(x1))<EPS_Newton))
{
cout<<"/nNewton method faild!";
return 0;
}
d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);
x0=x1;
cout<<"x("<<k<<")="<<x0<<endl;
}
while(fabs(d)>EPS_Newton&&fabs(f(x1))>ETA_Newton);
return 1;
}

double Newton(double (*f)(double),double (*f1)(double),double x0)
{
double x1,d;
int k=0;
do{
x1=x0-f(x0)/f1(x0);
if((k++>N_Newton)||(fabs(f1(x1))<EPS_Newton))
{
cout<<"/nNewton method faild!";
return 0;
}
d=(fabs(x1)<1?x1-x0:(x1-x0)/x1);
x0=x1;
cout<<"x("<<k<<")="<<x0<<endl;
}
while(fabs(d)>EPS_Newton&&fabs(f(x1))>ETA_Newton);
return 1;
}

float f_Newton(float x)
{
// return log(x)+x-2;
return x*x*x+x*x-3*x-3;
}

double f_Newton(double x)
{
// return log(x)+x-2;
return x*x*x+x*x-3*x-3;
}

float f1_Newton(float x)
{
// return 1/x+1;
return 3.0*x*x+2*x-3;
}

double f1_Newton(double x)
{
// return 1/x+1;
return 3.0*x*x+2*x-3;
}

void test_Newton()
{
float f_Newton(float);
float f1_Newton(float);
float x0,y0;
float Newton(float (*)(float),float (*)(float),float);

cout<<"Please input x0/n";
cin>>x0;
cout<<"x(0)="<<x0<<endl;
y0=Newton(f_Newton,f1_Newton,x0);
cout<<"/nThe root of the equation is x="<<(float)y0<<endl;
}

void test_Newton2()
{
double f_Newton(double);
double f1_Newton(double);
double x0,y0;
double Newton(double (*)(double),double (*)(double),double);

cout<<"Please input x0/n";
cin>>x0;
cout<<"x(0)="<<x0<<endl;
y0=Newton(f_Newton,f1_Newton,x0);
cout<<"/nThe root of the equation is x="<<(double)y0<<endl;
}
//////////////////////////////////////////////////////////////////////////
//牛顿值多项式
#define N_Difference 4
void Difference(float *x,float *y,int n)
{
float *f;
int k,i;
f =new float
;
for(k=1;k<=n;k++)
{
f[0]=y[k];
for(i=0;i<k;i++)
f[i+1]=(f[i]-y[i])/(x[k]-x[i]);
y[k]=f[k];
}
cout<<endl;
return ;
}

void Difference(double *x,double *y,int n)
{
double *f;
int k,i;
f =new double
;
for(k=1;k<=n;k++)
{
f[0]=y[k];
for(i=0;i<k;i++)
f[i+1]=(f[i]-y[i])/(x[k]-x[i]);
y[k]=f[k];
}
cout<<endl;
return ;
}

void test_Difference()
{
int i;
float varx =0.895,b;
float x[N_Difference+1] ={0.4,0.55,0.65,0.8,0.9};
float y[N_Difference+1] ={0.41075,0.57815,0.69675,0.88811,1.02652};
Difference(x,(float *)y,N_Difference);
b =y[N_Difference];
for (i=N_Difference-1;i>=0;i--)
{b = b*(varx-x[i])+y[i];
cout<<b<<endl;
}
cout<<"Nn("<<varx<<")="<<b<<endl;
}

void test_Difference2()
{
int i;
double varx =0.895,b;
double x[N_Difference+1] ={0.4,0.55,0.65,0.8,0.9};
double y[N_Difference+1] ={0.41075,0.57815,0.69675,0.88811,1.02652};
Difference(x,(double *)y,N_Difference);
b =y[N_Difference];
for (i=N_Difference-1;i>=0;i--)
{b = b*(varx-x[i])+y[i];
cout<<b<<endl;
}
cout<<"Nn("<<varx<<")="<<b<<endl;
}
//////////////////////////////////////////////////////////////////////////
//四阶阿当姆斯预测-校正公式
float *Runge_Kutta_Adams(float (*f)(float x,float y),float a,float b,float y0,int N)
{
float x=a,y=y0,K1,K2,K3,K4,*yy;
float h=(b-a)/N;
int i;
yy=new float[3];
for(i=1;i<=3;i++)
{
K1=f(x,y);
K2=f(x+h/2,y+h*K1/2);
K3=f(x+h/2,y+h*K2/2);
K4=f(x+h,y+h*K3);
y=y+h*(K1+2*K2+2*K3+K4)/6;
x=a+i*h;
*(yy+i-1)=y;
}
return yy;
}
double *Runge_Kutta_Adams(double (*f)(double x,double y),double a,double b,double y0,int N)
{
double x=a,y=y0,K1,K2,K3,K4,*yy;
double h=(b-a)/N;
int i;
yy=new double[3];
for(i=1;i<=3;i++)
{
K1=f(x,y);
K2=f(x+h/2,y+h*K1/2);
K3=f(x+h/2,y+h*K2/2);
K4=f(x+h,y+h*K3);
y=y+h*(K1+2*K2+2*K3+K4)/6;
x=a+i*h;
*(yy+i-1)=y;
}
return yy;
}
void Adams(float a,float b,int N,float (*f)(float x,float y),float y0)
{
int i;
float y1,y2,y,yp,yc,*yy,h,x;
cout<<"x[0]="<<a<<'/t'<<" y[0]="<<y0<<endl;
yy=Runge_Kutta_Adams(f,a,b,y0,N);
y1=yy[0];
y2=yy[1];
y=yy[2];
h=(b-a)/N;
for(i=1;i<=3;i++)
cout<<"x["<<i<<"]="<<a+i*h<<" y["<<i<<"]="<<*(yy+i-1)<<endl;
for(i=3;i<N;i++)
{
x=a+i*h;
yp=y+h*(55*f(x,y)-59*f(x-h,y2)+37*f(x-2*h,y1)-9*f(x-3*h,y0))/24.0;
yc=y+h*(9*f(x+h,yp)+19*f(x,y)-5*f(x-h,y2)+f(x-2*h,y1))/24.0;
cout<<"x["<<i+1<<"]="<<x+h<<" y["<<i+1<<"]="<<yc<<endl;
y0=y1;
y1=y2;
y2=y;
y=yc;
}
}
void Adams(double a,double b,int N,double (*f)(double x,double y),double y0)
{
int i;
double y1,y2,y,yp,yc,*yy,h,x;
cout<<"x[0]="<<a<<'/t'<<" y[0]="<<y0<<endl;
yy=Runge_Kutta_Adams(f,a,b,y0,N);
y1=yy[0];
y2=yy[1];
y=yy[2];
h=(b-a)/N;
for(i=1;i<=3;i++)
cout<<"x["<<i<<"]="<<a+i*h<<" y["<<i<<"]="<<*(yy+i-1)<<endl;
for(i=3;i<N;i++)
{
x=a+i*h;
yp=y+h*(55*f(x,y)-59*f(x-h,y2)+37*f(x-2*h,y1)-9*f(x-3*h,y0))/24.0;
yc=y+h*(9*f(x+h,yp)+19*f(x,y)-5*f(x-h,y2)+f(x-2*h,y1))/24.0;
cout<<"x["<<i+1<<"]="<<x+h<<" y["<<i+1<<"]="<<yc<<endl;
y0=y1;
y1=y2;
y2=y;
y=yc;
}
}

float f_Adams(float x,float y)
{
return -x*y*y;
}
double f_Adams(double x,double y)
{
return -x*y*y;
}
void test_Adams()
{
float a=0,b=5.0,y0=2.0;
int N=20;
Adams(a,b,N,f_Adams,y0);
}
void test_Adams2()
{
double a=0,b=5.0,y0=2.0;
int N=20;
Adams(a,b,N,f_Adams,y0);
}

//////////////////////////////////////////////////////////////////////////
//雅可比迭代法
//雅可比迭代法
double EPS_Jacobi =1e-6;
double MAX_Jacobi =100;
float *Jacobi(float a[3][4],int n)
{
float *x,*y,epsilon,s;
int i,j,k=0;

x =new float
; //here must be modifed
y =new float
; //here must be modifed
for( i=0;i<n;i++ )
x[i]=0;
while(1)
{
epsilon=0;
k++;
for( i=0;i<n;i++ )
{
s=0;
for( j=0;j<n;j++ )
{
if(j==i)
continue;
s+=a[i][j]*x[j];
}
y[i] =(a[i]
-s)/a[i][i];
epsilon +=fabs(y[i]-x[i]);
}
if(epsilon<EPS_Jacobi)
{
cout<<"迭代次数为："<<k<<endl;
return y;
}
if( k>=MAX_Jacobi )
{
cout<<"The method is disconvergent";
return y;
}
for( i=0;i<n;i++ )
x[i]=y[i];
}
return y;
}

double *Jacobi(double a[3][4],int n)
{
double *x,*y,epsilon,s;
int i,j,k=0;

x =new double
; //here must be modifed
y =new double
; //here must be modifed
for( i=0;i<n;i++ )
x[i]=0;
while(1)
{
epsilon=0;
k++;
for( i=0;i<n;i++ )
{
s=0;
for( j=0;j<n;j++ )
{
if(j==i)
continue;
s+=a[i][j]*x[j];
}
y[i] =(a[i]
-s)/a[i][i];
epsilon +=fabs(y[i]-x[i]);
}
if(epsilon<EPS_Jacobi)
{
cout<<"迭代次数为："<<k<<endl;
return y;
}
if( k>=MAX_Jacobi )
{
cout<<"The method is disconvergent";
return y;
}
for( i=0;i<n;i++ )
x[i]=y[i];
}
return y;
}

void test_Jacobi()
{
int i;
float a[3][4]={5.0,2.0,1.0,8.0,2.0,8.0,-3.0,21.0,1.0,-3.0,-6.0,1.0}; //here must be modifed
float *x;
x=new float[3]; //here must be modifed
x=Jacobi(a,3);
for( i=0;i<3;i++ )
cout<<"x["<<i<<"]="<<x[i]<<endl;
}

void test_Jacobi2()
{
int i;
double a[3][4]={5.0,2.0,1.0,8.0,2.0,8.0,-3.0,21.0,1.0,-3.0,-6.0,1.0}; //here must be modifed
double *x;
x=new double[3]; //here must be modifed
x=Jacobi(a,3);
for( i=0;i<3;i++ )
cout<<"x["<<i<<"]="<<x[i]<<endl;
}

//////////////////////////////////////////////////////////////////////////
//自适应梯形公式（变步长）
int n_AutoTrap;
float AutoTrap(float (*f)(float),float a,float b)
{
int i;
float x,s,h=b-a;
float t1,t2=h/2.0*(f(a)+f(b));
n_AutoTrap=1;
do{
s=0.0;
t1=t2;
for(i=0;i<=n_AutoTrap-1;i++)
{
x=a+i*h+h/2;
s+=f(x);
}
t2=(t1+s*h)/2.0;
n_AutoTrap*=2;
h/=2.0;
}while( fabs(t2-t1)>1e-6);
return t2;
}

double AutoTrap(double (*f)(double),double a,double b)
{
int i;
double x,s,h=b-a;
double t1,t2=h/2.0*(f(a)+f(b));
n_AutoTrap=1;
do{
s=0.0;
t1=t2;
for(i=0;i<=n_AutoTrap-1;i++)
{
x=a+i*h+h/2;
s+=f(x);
}
t2=(t1+s*h)/2.0;
n_AutoTrap*=2;
h/=2.0;
}while( fabs(t2-t1)>1e-6);
return t2;
}

float f_AutoTrap(float x)
{
return 1/(1+x*x);
}
//用变步长梯形法计算积分积分只要将函数定义为

/*float f_AutoTrap(float x)
{
if(x==0)
return 1;
else
return sin(x)/x;
}*/
double f_AutoTrap(double x)
{
return 1/(1+x*x);
}
//用变步长梯形法计算积分积分只要将函数定义为

/*double f_AutoTrap(double x)
{
if(x==0)
return 1;
else
return sin(x)/x;
}*/

void test_AutoTrap()
{
float s;
float AutoTrap(float (*)(float),float,float);
s=AutoTrap(f_AutoTrap,0.0,1.0);
cout<<"T("<<n_AutoTrap<<")="<<s<<endl;
}

void test_AutoTrap2()
{
double s;
double AutoTrap(double (*)(double),double,double);
s=AutoTrap(f_AutoTrap,0.0,1.0);
cout<<"T("<<n_AutoTrap<<")="<<s<<endl;
}
//////////////////////////////////////////////////////////////////////////
//最小二乘法
float power(int i,float v)
{
float a=1.0;
while(i--)
a*=v;
return a;
}
double power(int i,double v)
{
double a=1.0;
while(i--)
a*=v;
return a;
}

float *Approx(float *x,float *y,int m,int n)
{
float *c,*a;
int i,j,t;
float power(int ,float);
float *ColPivot(float *,int );
c= new float[(n+1)*(n+2)];
for(i=0;i<=n;i++)
{
for (j=0;j<=n;j++)
{
*(c+i*(n+2)+j)=0.0;
for(t=0;t<=m-1;t++)
*(c+i*(n+2)+j)+=power(i+j,x[t]);
}
*(c+i*(n+2)+n+1)=0.0;
for(j=0;j<=m-1;j++)
*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);
}
a=ColPivot((float *)c,n+1);
return a;
}

double *Approx(double *x,double *y,int m,int n)
{
double *c,*a;
int i,j,t;
double power(int ,double);
double *ColPivot(double *,int );
c= new double[(n+1)*(n+2)];
for(i=0;i<=n;i++)
{
for (j=0;j<=n;j++)
{
*(c+i*(n+2)+j)=0.0;
for(t=0;t<=m-1;t++)
*(c+i*(n+2)+j)+=power(i+j,x[t]);
}
*(c+i*(n+2)+n+1)=0.0;
for(j=0;j<=m-1;j++)
*(c+i*(n+2)+n+1)+=y[j]*power(i,x[j]);
}
a=ColPivot((double *)c,n+1);
return a;
}

void test_Approx()
{
int i;
float *a;
float x[16]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
float y[16]={4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86,
10.00,10.20,10.32,10.42,10.50,10.55,10.58,10.60};
float *Approx(float *,float *,int ,int );
a=Approx(x,y,16,2);
for(i=0;i<=2;i++)
cout<<"a["<<i<<"]="<<a[i]<<endl;
}

void test_Approx2()
{
int i;
double *a;
double x[16]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};
double y[16]={4.00,6.40,8.00,8.80,9.22,9.50,9.70,9.86,
10.00,10.20,10.32,10.42,10.50,10.55,10.58,10.60};
double *Approx(double *,double *,int ,int );
a=Approx(x,y,16,2);
for(i=0;i<=2;i++)
cout<<"a["<<i<<"]="<<a[i]<<endl;
}

//////////////////////////////////////////////////////////////////////////
} //

===========================================

我的E－mail：zzubegtostudy@163.com，没事别骚扰。

2006.10.08 发表

2006.10.09 修改