头文件:
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* nlfunc.h
*
* Nonlinear function.
*
* Zhang Ming, 2010-04, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef NLFUNC_H
#define NLFUNC_H
namespace splab
{
template <typename Type>
class NLFunc
{
public:
/**
* Initialize the parameters
*/
NLFunc( Type aa, Type bb, Type cc ) : a(aa), b(bb), c(cc)
{ }
/**
* Compute the value of objective function at point x.
*/
Type operator()( Type &x )
{
return ( a*x*x + b*x + c );
}
/**
* Compute the gradient of objective function at point x.
*/
Type grad( Type &x )
{
return ( 2*a*x + b );
}
private:
// parameters
Type a,
b,
c;
};
// class NLFunc
}
// namespace splab
#endif
// NLFUNC_H
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* nleroot.h
*
* Solution for nonlinear equation.
*
* This file includes three routines for finding root of nonlinear equation.
* The bisection method don't need compute the function's gradient, but it
* has a slow convergence rate, linear convergence. Newton method has a
* quadratic convergence, however, it has to compute the gradient of function.
* The secant method don't need compute the gradient, in adition, it has a
* superlinear convergence rate( the order is 1.618 ), so it's a practical
* method in many cases.
*
* Zhang Ming, 2010-04, Xi'an Jiaotong University.
*****************************************************************************/
#ifndef NLEROOT_H
#define NLEROOT_H
#include <cmath>
#include <constants.h>
#include <nlfunc.h>
using namespace std;
namespace splab
{
template<typename Type> Type bisection( NLFunc<Type> &f, Type a, Type b,
Type tol=Type(1.0e-6) );
template<typename Type> Type newton( NLFunc<Type> &f, Type x0,
const Type tol=Type(1.0e-6),
int maxItr=MAXTERM );
template<typename Type> Type secant( NLFunc<Type> &f, Type x1, Type x2,
const Type tol=Type(1.0e-6),
int maxItr=MAXTERM );
#include <nleroot-impl.h>
}
// namespace splab
#endif
// NLEROOT_H
实现文件:
/*
* Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
*
* This program is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License as published by the
* Free Software Foundation, either version 2 or any later version.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
*
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* This program is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
* more details. A copy of the GNU General Public License is available at:
* http://www.fsf.org/licensing/licenses
*/
/*****************************************************************************
* nleroot-impl.h
*
* Implementation for nonlinear equation rooting.
*
* Zhang Ming, 2010-04, Xi'an Jiaotong University.
*****************************************************************************/
/**
* Bisection method for finding function root.
*/
template <typename Type>
Type bisection( NLFunc<Type> &f, Type a, Type b, Type tol )
{
Type c = 0,
fc = 0,
fa = f(a),
fb = f(b);
if( a > b )
{
c = a;
a = b;
b = c;
}
if( fa*fb > 0 )
{
cerr << " Invalid interval!" << endl;
return Type(0);
}
int maxItr = (int) ceil( (log(b-a)-log(tol)) / log(2.0) );
for( int i=0; i<maxItr; ++i )
{
c = ( a+ b) / 2;
fc = f(c);
if( abs(fc) < EPS )
return c;
if( fb*fc > 0 )
{
b = c;
fb = fc;
}
else
{
a = c;
fa = fc;
}
if( (b-a) < tol )
return (b+a)/2;
}
cout << "No solution for the specified tolerance!" << endl;
return c;
}
/**
* Newton method for finding function root.
*/
template <typename Type>
Type newton( NLFunc<Type> &f, Type x0, Type tol, int maxItr )
{
Type xkOld = x0,
xkNew,
fkGrad;
for( int i=0; i<maxItr; ++i )
{
fkGrad = f.grad(xkOld);
while( abs(fkGrad) < EPS )
{
xkOld *= 1.2;
fkGrad = f.grad(xkOld);
}
xkNew = xkOld - f(xkOld)/fkGrad;
if( abs(f(xkNew)) < tol )
return xkNew;
if( abs(xkOld-xkNew) <tol )
return xkNew;
xkOld = xkNew;
}
cout << "No solution for the specified tolerance!" << endl;
return xkNew;
}
/**
* Secant method for finding function root.
*/
template <typename Type>
Type secant( NLFunc<Type> &f, Type x1, Type x2, Type tol, int maxItr )
{
Type x_k,
x_k1 = x1,
x_k2 = x2,
f_k1 = f(x_k1),
f_k2;
for( int i=0; i<maxItr; ++i )
{
f_k2 = f(x_k2);
if( abs(f_k2-f_k1) < EPS )
{
x_k2 = (x_k1+x_k2) / 2;
f_k2 = f(x_k2);
}
x_k = x_k2 - (x_k2-x_k1)*f_k2/(f_k2-f_k1);
if( abs(f(x_k)) < tol )
return x_k;
if( abs(x_k2-x_k) < tol )
return x_k;
f_k1 = f_k2;
x_k1 = x_k2;
x_k2 = x_k;
}
cout << "No solution for the specified tolerance!" << endl;
return x_k;
}
测试代码:
/*****************************************************************************
* nle_test.cpp
*
* Rooting of nonlinear equation testing.
*
* Zhang Ming, 2010-04, Xi'an Jiaotong University.
*****************************************************************************/
#include <iostream>
#include <iomanip>
#include <nleroot.h>
using namespace std;
using namespace splab;
typedef double Type;
int main()
{
Type x01, x02, x03;
NLFunc<Type> f( 1.0, -3.0, 1.0 );
cout << setiosflags(ios::fixed) << setprecision(8);
x01 = bisection( f, 0.0, 2.0, 1.0e-6 );
cout << "Bisection method : " << x01 << endl <<endl;
x02 = newton( f, 0.0, 1.0e-6, 100 );
cout << "Newton method : " << x02 << endl <<endl;
x03 = secant( f, 1.0, 3.0, 1.0e-6, 100 );
cout << "Secant method : " << x03 << endl <<endl;
return 0;
}
运行结果:
Bisection method : 0.38196611
Newton method : 0.38196601
Secant method : 2.61803406
Process returned 0 (0x0) execution time : 0.062 s
Press any key to continue.