Actual code:
Code:
//Newton-Raphson Method
//http://en.wikipedia.org/wiki/Newton's_method
#include <cstdlib>
#include <iostream>
#include <math.h>
using namespace std;
//------------------------------------------------------------------------------
class point
{
public:
double x1, x2;
double length()
{
return sqrt(pow(x1,2)+pow(x2,2)); //sqrt(x1^2+x2^2)
}
//I write that overloaded operator is a friend of this class.
friend ostream & operator<<(ostream & theStream, point & thePoint);
};
//overloading operators point*double, double*point, point+point, point-point
point operator*(point x, double a)
{
point result;
result.x1 = x.x1 * a;
result.x2 = x.x2 * a;
return result;
}
point operator*(double a, point x)
{
point result;
result.x1 = x.x1 * a;
result.x2 = x.x2 * a;
return result;
}
point operator+(point & thePoint1, point & thePoint2)
{
point sum;
sum.x1 = thePoint1.x1 + thePoint2.x1;
sum.x2 = thePoint1.x2 + thePoint2.x2;
return sum;
}
point operator-(point & thePoint1, point & thePoint2)
{
point sum;
sum.x1 = thePoint1.x1 - thePoint2.x1;
sum.x2 = thePoint1.x2 - thePoint2.x2;
return sum;
}
//overloaded operator <<
ostream & operator<<(ostream & theStream, point & thePoint)
{
theStream<<"["<<thePoint.x1<<", "<<thePoint.x2<<"]";
return theStream;
}
//------------------------------------------------------------------------------
class hessian
{
public: //I've got a problem with determinant (in main function there is
//a problem with compiling if I use 'private' here!!!!!!!!
double x00, x01, x10, x11;
//|x00 x01|
//|x10 x11|
double determinant;
public:
friend point operator*(hessian theHessian, point thePoint);
friend hessian operator*(hessian theHes1, hessian theHes2);
void det();
void invertHessian();
void hessianAtPoint(point thePoint);
void hessianAtPoint();
friend ostream&operator<<(ostream & theStream, hessian & theHessian);
};
//internal functions - determinant, invert hessian and hessian at point
void hessian::det()
{ determinant = x00 * x11 - x10 * x01; }
void hessian::invertHessian()
{
//inverted matrix: A^(-1) = A^D / det(A)
//A = [(a, b), (c, d)], A^D = [(d, -b), (-c, a)]
double _x00, _x01, _x10, _x11;
_x00 = x11/determinant;
_x01 = -x01/determinant;
_x10 = x10/determinant;
_x11 = x00/determinant;
x00 = _x00; x01 = _x01; x10 = _x10; x11 = _x11;
}
void hessian::hessianAtPoint()
{
x00 = x11 = 2;
x01 = x10 = 0;
}
//overloading operators << and hessian*hessian
ostream&operator<<(ostream & theStream, hessian & theHessian)
{
theStream<<"|"<<theHessian.x00<<" "<<theHessian.x01<<"|]n";
theStream<<"|"<<theHessian.x10<<" "<<theHessian.x11<<"|]n";
return theStream;
}
hessian operator*(hessian theHes1, hessian theHes2)
{
hessian result;
//simple multiplication matrix by matrix
//[a b][e f] = [ae+bg af+bh]
//[c d][g h] [ce+dg cf+dh]
result.x00 = theHes1.x00 * theHes2.x00 + theHes1.x01 * theHes2.x10;
result.x01 = theHes1.x00 * theHes2.x01 + theHes1.x01 * theHes2.x11;
result.x10 = theHes1.x10 * theHes2.x00 + theHes1.x11 * theHes2.x10;
result.x11 = theHes1.x10 * theHes2.x01 + theHes1.x11 * theHes2.x11;
return result;
}
//------------------------------------------------------------------------------
class gradient
{
public: //???????????????????
point itsPoint;
public:
void gradientAtPoint(point thePoint);
friend ostream&operator<<(ostream & theStream, gradient & theGradient);
};
void gradient::gradientAtPoint(point thePoint)
{
itsPoint.x1 = 2 * thePoint.x1 - 1;
itsPoint.x2 = 2 * thePoint.x2 - 1;
}
ostream&operator<<(ostream & theStream, gradient & theGradient)
{
theStream<<theGradient.itsPoint;
return theStream;
}
//------------------------------------------------------------------------------
//overloading other operator: hessian*point
point operator*(hessian theHessian, point thePoint)
{
point result;
//[x00 x01] * [x1 x2] = [x00*x1 + x01*x2]
//[x10 x11] [x10*x1 + x11*x2]
result.x1 = theHessian.x00 * thePoint.x1 + theHessian.x01 * thePoint.x2;
result.x2 = theHessian.x10 * thePoint.x1 + theHessian.x11 * thePoint.x2;
return result;
}
//main function
double myFunction(point thePoint)
{ return (pow(thePoint.x1, 2)+pow(thePoint.x2, 2)); };
//------------------------------------------------------------------------------
int main()
{
//declaring variables
gradient myGradient;
hessian myHessian;
point myPoint;
point tempPoint;
double accuracy;
int n = 0;
//stating input variables
cout<<"Write the input value x1 = ";
cin>>myPoint.x1;
cout<<"Write the input value x2 = ";
cin>>myPoint.x2;
cout<<"Write the accuracy = ";
cin>>accuracy;
//main algorithm
myGradient.gradientAtPoint(myPoint);
while (myGradient.itsPoint.length() > accuracy)
{
n++;
myGradient.gradientAtPoint(myPoint);
//myHessian.hessianAtPoint(myPoint);
myHessian.hessianAtPoint();
myHessian.det();
myHessian.invertHessian();
//myPoint = myPoint - (myHessian * (myGradient.itsPoint)); //[1]
tempPoint = myHessian * myGradient.itsPoint; //[2]
myPoint = myPoint - tempPoint; //[3]
//Here I changed line [1] into lines [2] and [3].
//I think both [1] and [2][3] are the same but [1] doesn't compile.
}
//writing output
cout<<"\n[x1, x2] = ["<<myPoint.x1<<", "<<myPoint.x2<<"]";
cout<<"\nminimum value = "<<myFunction(myPoint);
cout<<"\naccuracy = "<<accuracy;
cout<<"\nnumber of iterations = "<<n<<"\n";
system("PAUSE");
return EXIT_SUCCESS;
}
//And by the way - I've got other question. If I'd like the user to be able to
//change the function, is there any way to allow this application to read the
//fuction (and gradient, hessian) from the file and put inside those lines which
//are responsible for calculations?
1. I changed
Code:
int main(int argc, char *argv[])
into
2. I decided to overload hessianAtPoint. I don't use hessianAtPoint(point thePoint) at all but I use hessianAtPoint.
Code:
void hessianAtPoint(point thePoint);
void hessianAtPoint();
Code:
void hessian::hessianAtPoint()
{
x00 = x11 = 2;
x01 = x10 = 0;
}
Code:
//myHessian.hessianAtPoint(myPoint);
myHessian.hessianAtPoint();
You asked why I take myPoint into that function if I don't use it at all. It is because this program should calculate something. As "input" you've got function - here it is
Code:
//main function
double myFunction(point thePoint)
{ return (pow(thePoint.x1, 2)+pow(thePoint.x2, 2)); };
Then, with using calculations on the sheet of paper, I found hessian at point (matrix of partial derivatives) and gradient at point (some other partial derivatives). If the function is easy, then hessian and gradient may contain only numbers, like here:
Code:
x00 = x11 = 2;
x01 = x10 = 0;
But if function is more complicated, there may be other function dependent on thePoint, for example:
Code:
void HessianAtPoint(point c)
{
//x00=x11=2;
//x01=x10=0;
x00=(((-2*((-7.9)*exp(-7.9*c.x1)+0.1*exp(0.1*c.x1))*pow(c.x2,2))/pow((exp(-7.9*c.x1)+exp(0.1*c.x1)),2))+((-2*((-7.9)*exp(-7.9*c.x1)+0.1*exp(0.1*c.x1))*(-1+c.x1*pow(c.x2,2)))/pow((exp(-7.9*c.x1)+exp(0.1*c.x1)),3))-((((62.41)*exp(-7.9*c.x1)+0.01*exp(0.1*c.x1))*(-1+c.x1*pow(c.x2,2)))/pow((exp(-7.9*c.x1)+exp(0.1*c.x1)),2)));
x01=x10=((2*c.x2*(exp((-7.9)*c.x1)+exp(0.1*c.x1))-2*c.x1*c.x2*((-7.9)*exp((-7.9)*c.x1)+0.1*exp(0.1*c.x1)))/(pow((exp(-7.9*c.x1)+exp(0.1*c.x1)),2)));
x11=((2*c.x1)/(exp(-7.9*c.x1)+exp(0.1*c.x1)));
}
I don't know how to create "dynamic" which take formula of function, formulas of hessian at point and formulas of gradient at point. This is why I decided to write it inside the code and not to take it as a parameter from the user. So if I want to use this program for other function, I need to change some lines inside the code and recompile it.
3. I added return:
Code:
ostream&operator<<(ostream & theStream, hessian & theHessian)
{
theStream<<"|"<<theHessian.x00<<" "<<theHessian.x01<<"|]n";
theStream<<"|"<<theHessian.x10<<" "<<theHessian.x11<<"|]n";
return theStream;
}
4. I don't understand this error
Error 5 error C4716: 'operator*' : must return a value g:\w00t\visual studio 2008\projects\temp\temp2.cpp 126
I have got:
Code:
hessian operator*(hessian theHes1, hessian theHes2)
{
hessian result;
//simple multiplication matrix by matrix
//[a b][e f] = [ae+bg af+bh]
//[c d][g h] [ce+dg cf+dh]
result.x00 = theHes1.x00 * theHes2.x00 + theHes1.x01 * theHes2.x10;
result.x01 = theHes1.x00 * theHes2.x01 + theHes1.x01 * theHes2.x11;
result.x10 = theHes1.x10 * theHes2.x00 + theHes1.x11 * theHes2.x10;
result.x11 = theHes1.x10 * theHes2.x01 + theHes1.x11 * theHes2.x11;
return result;
}
Result is hessian, hessian has got four internal variables:
Code:
class hessian
{
public: //I've got a problem with determinant (in main function there is
//a problem with compiling if I use 'private' here!!!!!!!!
double x00, x01, x10, x11;
So I'm just writing resulting x00, x01, x10 and x11 into object which is hessian. To tell the truth I don't remember where I used (or whether I used) this overloaded operator to multiply hessian by hessian so probably I don't need to overload this operator hessian*hessian and I can simply remove this part of code.
5. I added
return result;
into
Code:
//overloading other operator: hessian*point
point operator*(hessian theHessian, point thePoint)
{
point result;
//[x00 x01] * [x1 x2] = [x00*x1 + x01*x2]
//[x10 x11] [x10*x1 + x11*x2]
result.x1 = theHessian.x00 * thePoint.x1 + theHessian.x01 * thePoint.x2;
result.x2 = theHessian.x10 * thePoint.x1 + theHessian.x11 * thePoint.x2;
return result;
}
Regards
Originally Posted by zdzislavv
. I guess there are some other things which may be improved in this code (for example proper way of choosing which variables are private and which are public in classes). Are these classes constructed properly or should I change something?
