数值分析源码和库

最新推荐文章于 2022-11-28 23:13:26 发布
最新推荐文章于 2022-11-28 23:13:26 发布
阅读量2.8k 收藏 5
点赞数 1
CC 4.0 BY-SA版权
分类专栏： 源码模块/Source Code 文章标签： float c 算法 ini delete input
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/begtostudy/article/details/1326534

编前语：

这次放出的是只有源码，大部分是拼凑来的，我都做了修改（郁闷那么多错还发到网上共享），多为通用算法，专业级的算法有，乘我还是学生我会想办法为大家弄的。以后会放出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[n];
 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[n];
 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[n];
 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[n];
 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[n];
 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[n];
 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[n]; 
 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[n]; 
 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[n]; //here must be modifed
 y =new float[n]; //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][n]-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[n]; //here must be modifed
 y =new double[n]; //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][n]-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 修改 

数值分析算法源代码
08-15
数值分析一些常用算法的MATLAB源代码，包括二分法、newton法等15个常用算法代码。
数值分析方法库
11-04
这是一个关于数值分析的一个方法库，主要是程序的编写
用于数值计算的开源库 GNU Scientific Library
做好程序员
09-08 3330
一直很好奇使用C/C++语言，该如何进行复杂的数值(数学)计算：例如解方程组，C的标准数学库就做不到了。 最近刚了解到, 这些问题已经有很好用的开源库了，它就是： GSL - GNU Scientific Libraryhttp://www.gnu.org/softwar
刘寅立matlab数值计算案例分析 pdf,《MATLAB数值计算案例分析》源代码
weixin_32771227的博客
03-17 749
function [x, y, exitflag] = GoldenSectionFun(fun, a, b, eps, maxstep)% 黄金分割法优化函数% 输入参数：% fun为优化函数% a为取值区间左侧值% b为取值区间右侧值% eps为精度约束% maxstep为最大迭代次数约束% 输出参数：% x为最优解% y为最小值% exitflag为退出标记% 输入参数检查i...
基于Matlab实现常用数值分析程序（源码）.rar
06-15
在数学与工程计算领域，Matlab是一款不可或缺的工具，它以其强大的矩阵运算能力和丰富的内置函数库，成为科研人员和工程师进行数值分析的首选平台。本资源集合提供了基于Matlab实现的常用数值分析程序源码，旨在帮助...
基于Python的数值分析算法实现设计源码
09-27
此外，该源码中还包含了与数值分析相关的科学计算库文件，比如numpy和matplotlib。numpy是Python中用于进行高效数值计算的基础库，它提供了多维数组对象以及一系列处理这些数组的函数。matplotlib则是一个用于创建...
作业_数值分析matlab程序_源码
10-04
在数值分析领域，MATLAB是一种广泛使用的工具，它提供了丰富的函数库和强大的计算能力，使得复杂的数值计算变得简单。这份“作业_数值分析matlab程序_源码”压缩包文件显然是针对学习数值分析的学生或者研究人员准备...
基于Python和MATLAB的电容测量电路数值仿真源码设计与对比分析
最新发布
10-02
本文主要探讨了基于Python和MATLAB的电容测量电路数值仿真源码的设计与对比分析。项目中包括了MATLAB和Python两种语言编写的仿真程序，旨在通过这两种不同平台的对比，验证Python在实现电容测量电路仿真中的可行性、...
C#数值分析编程源码
12-12
这个名为"C#数值分析编程源码"的压缩包很显然包含了使用C#实现的各种数值分析算法的源代码。数值分析是应用数学的一个分支，它涉及到解决实际问题中的数值计算方法，例如求解微积分方程、优化问题、傅里叶变换等。 ...
数值分析——牛顿迭代法c++源代码
04-17
数值分析——牛顿迭代法c++源代码.程序十分简洁明了，相信会对你有帮助
数值分析（经典算法）
10-13
学习数值分析的人，里面的算法都是经常必备的常用经典算法。
数值分析的迭代mathematic代码
04-16
Fixed-Point Iteration Newton's Method
数值计算方法 雅克比迭代法实验报告及源码
05-02
数值计算方法课程中的一个上机实验题目 雅克比迭代，一份很详细的实验报告，里面包括了源码
科学计算与数值分析方面的算法库
love code, love life.
11-28 879
科学计算和数值分析是一个研究领域，旨在提供借助计算机解决不同科学领域的大规模问题的方法。典型的问题是常微分方程和偏微分方程（ODE，PDE），它们的离散化，以及由此产生的线性代数问题的解。这里列出了常见的资源，包括基本的线性代数库，多目的工具箱，有限元，网格化，数据格式，稀疏线性求解，可视化，其他库和工具，以及社区。这些库涵盖了多种语言，如C/C++，python等。
数模算法-数值分析算法
Dompink的博客
09-03 1824
数值分析（numerical analysis），是数学的一个分支，主要研究连续数学（区别于离散数学）问题的 算法。 如果在比赛中采用高级语言进行编程的话，那一些数值分析中常用的算法比 如方程组求解、矩阵运算、 函数积分等算法就需要额外编写库函数进行调用。 这类算法是针对高级语言而专门设的，如果你用的是 MATLAB 、 Mathematica ，大可不必准备， 因为像数值分析中有很多函数
常用数值计算库
~。~|||
10-28 1345
 program 常用数值计算库Intel Math Kernel Library1.基本线形代数运算(BLAS) 向量与向量、向量与矩阵、矩阵与矩阵的运算2.稀疏线形代数运算3.快速傅立叶变换(单精度/双精度)(fftw)4.LAPACK(求解线形方程组、最小方差、特征值、Sylvester方程等)5.向量数学库(VML)6.向量统计学库(VSL)7.高级离散傅立叶变换IMSL软件名称 IMSL
数值分析多种算法C语言代码－推荐
worldy的专栏
12-10 3368
数值分析多种算法C语言代码－推荐 分类： DSP数字信号处理 2013-10-15 22:20 157人阅读 评论(0) 收藏 举报 1、离散傅立叶变换与反变换 [cpp] view plaincopyprint? /****************************************************************
MATLAB数值分析源码章节整合指南
每个章节都对应着数值分析中的一个重要主题，并且在MATLAB中都有相对应的算法实现和函数库支持。 通过深入研究这些源码，可以了解到MATLAB在不同数值分析问题中算法实现的细节，掌握其适用场景及优缺点。这对于提高...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值