Минимизация Методом Ньютона уравнения с двумя переменными

#include "stdafx.h"
#include <math.h>
#include <conio.h>
#include <stdio.h>
#include <iostream>
 
 
double f(double x)
{
    double f;
    f = exp(-x)*log(x);
    return f;
}
double df(double x)
{
    double f;
    f = -exp(-x)*log(x) + (exp(-x) / x);
    return f;
}
 
void main()
{
    double e;
    std::cout << "Vvedite e" << std::endl;
    std::cin >> e;  
    double y1, y2, z1, z2, s, z, y, x, a = 1, b = 3;
    int N;
    //N=0; //Schetchik iteracii
    y1 = f(a);
    y2 = f(b);
    z1 = df(a);
    z2 = df(b);
    N = 4;
    do
    {
        s = ((z2*b - z1*a) - (y2 - y1)) / (z2 - z1);
        y = f(s);
        z = df(s);
        N += 2;
        if (z == 0) { x = s; break; }
        if (z>0) { b = s; y2 = y; z2 = z; }
        else { a = s; y1 = y; z1 = z; }
    } while ((b - a)/2>e);
    if (z != 0) { x = (a + b) / 2; y = f(x); }
    printf("Metod Newtona-Rafsona \n \n");
    printf("x = %.5f \n", x);
    printf("y = %.5f \n \n", y);
    printf("Chislo iteracii: %i", N);
}

#include <iostream>
#include <vector>
#include <ctime>
#include <algorithm>
#include <iterator>
#include <cmath>
#include "simplex.h"
using namespace std;
 
template<class Con>
void printcon(const Con& c){
  std::cout.precision(12);
  copy( c.begin(), 
    c.end(), 
    ostream_iterator<typename Con::value_type>(cout, "  ") );
  cout<<endl;
}
 
class Rosenbrock{
public:
  double  operator()(vector<double> x){
    return 100*pow(x[1]-pow(x[0],2),2)+pow(1-x[0],2);
  }
};
 
double rosenbrock(vector<double> x){
  return 100*pow(x[1]-pow(x[0],2),2)+pow(1-x[0],2);
}
 
double polynomial_fun(vector<double> x){
  return pow( pow(x[0],2) + pow(x[1],2), 2 ) - pow(x[0]-3*x[1],2);
}
 
int main(){
 
  //initial guessing value for Rosenbrock function
  vector<double> init;
  init.push_back(1.23);
  init.push_back(10.96);
 
  //initial trial simplex
  double a0[]={-1.5, -1}, a1[]={-2, -2}, a2[]={ 2.5, 1.5};
  vector<vector<double> > simplex;
  simplex.push_back( vector<double>(a0,a0+2) );
  simplex.push_back( vector<double>(a1,a1+2) );
  simplex.push_back( vector<double>(a2,a2+2) );
 
  //optimizating...
  //printcon is a function printing container values
 
  using BT::Simplex;
  cout<<"Rosenbrock function achieves minimum at:"<<endl;
  printcon( Simplex(Rosenbrock(), init) );  cout<<endl;
  cout<<"Rosenbrock function achieves minimum at:"<<endl;
  printcon( Simplex(rosenbrock, init, 1e-7, simplex, 1E5) ); cout<<endl;
  cout<<"Polynomial function achieves minimum at:"<<endl;
  printcon( Simplex(polynomial_fun, init, 1e-7, simplex) );  cout<<endl;
}

#include <vector>
#include <limits>
#include <algorithm>
#include <functional>
#include <iostream>
 
namespace BT{
 
  template<class D, class OP>
  std::vector<D> Simplex(OP f,                   //target function
             std::vector<D> init,    //initial guess of the parameters
             D tol=1E8*std::numeric_limits<D>::epsilon(), //termination criteria
             std::vector<std::vector<D> > x =  std::vector<std::vector<D> >(),
             //x: The Simplex
             int iterations=1E5){    //iteration step number
 
    int N=init.size();                         //space dimension
    const double a=1.0, b=1.0, g=0.5, h=0.5;   //coefficients
                                               //a: reflection  -> xr  
                                               //b: expansion   -> xe 
                                               //g: contraction -> xc
                                               //h: full contraction to x1
    std::vector<D> xcentroid_old(N,0);   //simplex center * (N+1)
    std::vector<D> xcentroid_new(N,0);   //simplex center * (N+1)
    std::vector<D> vf(N+1,0);            //f evaluated at simplex vertexes       
    int x1=0, xn=0, xnp1=0;         //x1:   f(x1) = min { f(x1), f(x2)...f(x_{n+1} }
                                    //xnp1: f(xnp1) = max { f(x1), f(x2)...f(x_{n+1} }
                                    //xn:   f(xn)<f(xnp1) && f(xn)> all other f(x_i)
    int cnt=0; //iteration step number
 
    if(x.size()==0) //if no initial simplex is specified
      { //construct the trial simplex
    //based upon the initial guess parameters
      std::vector<D> del( init );
      std::transform(del.begin(), del.end(), del.begin(), 
             std::bind2nd( std::divides<D>() , 20) );//'20' is picked 
                                                             //assuming initial trail close to true
      
      for(int i=0; i<N; ++i){
    std::vector<D> tmp( init );
    tmp[i] +=  del[i];
    x.push_back( tmp );
      }
      x.push_back(init);//x.size()=N+1, x[i].size()=N
      
      //xcentriod
      std::transform(init.begin(), init.end(), 
        xcentroid_old.begin(), std::bind2nd(std::multiplies<D>(), N+1) );
      }//constructing the simplex finished
    
    //optimization begins
    for(cnt=0; cnt<iterations; ++cnt){
 
      for(int i=0;i<N+1;++i){
    vf[i]= f(x[i]);
      }
      
      x1=0; xn=0; xnp1=0;//find index of max, second max, min of vf.
      
      for(int i=0;i<vf.size();++i){
    if(vf[i]<vf[x1]){
      x1=i;
    }
    if(vf[i]>vf[xnp1]){
      xnp1=i;
    }
      }
      
      xn=x1;
      
      for(int i=0; i<vf.size();++i){ 
    if( vf[i]<vf[xnp1] && vf[i]>vf[xn] )
      xn=i;
      }
      //x1, xn, xnp1 are found
 
      std::vector<D> xg(N, 0);//xg: centroid of the N best vertexes
      for(int i=0; i<x.size(); ++i){
    if(i!=xnp1)
      std::transform(xg.begin(), xg.end(), x[i].begin(), xg.begin(), std::plus<D>() );
      }
      std::transform(xg.begin(), xg.end(), 
        x[xnp1].begin(), xcentroid_new.begin(), std::plus<D>());
      std::transform(xg.begin(), xg.end(), xg.begin(), 
        std::bind2nd(std::divides<D>(), N) );
      //xg found, xcentroid_new updated
 
      //termination condition
      D diff=0;          //calculate the difference of the simplex centers
                         //see if the difference is less than the termination criteria
      for(int i=0; i<N; ++i)     
    diff += fabs(xcentroid_old[i]-xcentroid_new[i]);
 
      if (diff/N < tol) break;              //terminate the optimizer
      else xcentroid_old.swap(xcentroid_new); //update simplex center
      
      //reflection:
      std::vector<D> xr(N,0); 
      for( int i=0; i<N; ++i)
    xr[i]=xg[i]+a*(xg[i]-x[xnp1][i]);
      //reflection, xr found
      
      D fxr=f(xr);//record function at xr
      
      if(vf[x1]<=fxr && fxr<=vf[xn])
    std::copy(xr.begin(), xr.end(), x[xnp1].begin() );
      
      //expansion:
      else if(fxr<vf[x1]){
    std::vector<D> xe(N,0);
    for( int i=0; i<N; ++i)
      xe[i]=xr[i]+b*(xr[i]-xg[i]);
    if( f(xe) < fxr )
      std::copy(xe.begin(), xe.end(), x[xnp1].begin() );
    else
      std::copy(xr.begin(), xr.end(), x[xnp1].begin() );
      }//expansion finished,  xe is not used outside the scope
      
      //contraction:
      else if( fxr > vf[xn] ){
    std::vector<D> xc(N,0);
    for( int i=0; i<N; ++i)
      xc[i]=xg[i]+g*(x[xnp1][i]-xg[i]);
    if( f(xc) < vf[xnp1] )
      std::copy(xc.begin(), xc.end(), x[xnp1].begin() );
    else{
      
      for( int i=0; i<x.size(); ++i ){
        if( i!=x1 ){ 
          for(int j=0; j<N; ++j) 
        x[i][j] = x[x1][j] + h * ( x[i][j]-x[x1][j] );
        }
      }
      
    }
      }//contraction finished, xc is not used outside the scope
 
    }//optimization is finished
 
    if(cnt==iterations){//max number of iteration achieves before tol is satisfied
      std::cout<<"Iteration limit achieves, result may not be optimal"<<std::endl;
    }
    return x[x1];
  }
}