/* Duel numbers consists of a real part x and an infinitesimal part y e. Where e is defined such that e^2 = 0 but e != 0 Hence let w = x + e*y and q = a + e*b Addition: w+q = (x+e*y)+(a+e*b) = x+e*y+a+e*b = (x+a) + (y+b)*e Substraction: w-q = (x+e*y)-(a+e*b) = x+e*y-a-e*b = (x-a) + (y-b)*e Product: w*q = (x+e*y)*(a+e*b) = x*a + e*x*b + a*e*y + e*y*e*b = x*a + (x*b + e*y)*e Division: w x+e*y x+e*y a-e*b x*a - x*b*e + a*y*e - y*b*e^2 --- = -------- = -------- * -------- = -------------------------------- q a+e*b a+e*b a-e*b a*a -a*b*e + a*b*e - b^2 e^2 x*a + (y*a-b*x) e = -------------------- a^2 NOTE: Taylor expansion for real function: f: R -> R : f(x+dx) = f(x) + f'(x)dx + f''(x)dx^2/2 + f'''(x)dx^3/6 + ... Taylor expansion for function of Dual number: f: D -> D : use dx=e*y f(x+e y) = f(x) + f'(x) e y + f''(x) (ey)^2/2 + f'''(x)(e y)^3/6 + ... = f(x) + f'(x) e y because e^2 = 0 Hence f(x+1e) = f(x) + f'(x) e We can prove all rules (f+g)' = f'(x)+g'(x) (fg)' = f'(x)g(x)+f(x)g'(x) (f/g)' = .. by Taylor expansion f(x)g(x) and collecting terms of different order. Lets prove (fg)' by Taylor expansion: f(x)g(x) = (f(x) + f'(x)dx + f''(x)dx^2/2)*( g(x) + g'(x)dx + g''(x)dx^2/2 ) = f(x)g(x) + ( f'(x)g(x)+f(x)g'(x) )dx + ( ... )dx^2 + .. Hence for dual numbers: f(x+ey)g(x+ey) = f(x)g(x) + ( f'(x)g(x)+f(x)g'(x) )e no higher order terms CODE details: Because all methods are just a few lines of code, I put everything inside the class definition, instead of defining a DualNumber.cpp file. NOTE: I only define some of the mathematical operators and functions. Hence you will get compile errors if you try to use operators / functions that I have not defined. */ #include #include using namespace std; class DualNumber { private: // internal variables that are not accessible outside the class. double real_value; double infitesimal_value; public: // These are contructors for declaring a dual number DualNumber(double r, double i) : real_value(r), infitesimal_value(i) { ; } DualNumber(double r) : real_value(r), infitesimal_value(0.0) { ; } DualNumber() : real_value(1.0), infitesimal_value(0.0) { ; } // here are get methods that returns the value of the internal variables. double real() const { return real_value; } double infi() const { return infitesimal_value; } DualNumber operator ()(const int& r) // Lift an integer to the corresponding dual number. e.g. used for DualNumber(42) { return DualNumber( double(r), 0 ); } DualNumber operator()(const double &r) // Lift a real number to the corresponding dual number. e.g. used for DualNumber(42.0) { return DualNumber( r, 0 ); } }; // Class ends here /* We need to define the meaning of mathematical operations between pairs of dual numbers. Fortunately, mathematical operators are functions in C++ that we can define ourselves. */ inline DualNumber operator+(const DualNumber &w1, const DualNumber &w2) // w+q { double r=w1.real() + w2.real(); double i=w1.infi() + w2.infi(); return DualNumber( r, i ); } inline DualNumber operator-(const DualNumber &w1, const DualNumber &w2) // w-q { double r=w1.real() - w2.real(); double i=w1.infi() - w2.infi(); return DualNumber( r, i ); } inline DualNumber operator*(const DualNumber &w1, const DualNumber &w2) // w*q { double r=w1.real()*w2.real(); double i=w1.real()*w2.infi()+w1.infi()*w2.real(); return DualNumber( r, i ); } inline DualNumber operator*(double& a, const DualNumber &w) // a*w { double r=a*w.real(); double i=a*w.infi(); return DualNumber( r, i ); } inline DualNumber operator*(const DualNumber &w, double& a) // w*a { double r=a*w.real(); double i=a*w.infi(); return DualNumber( r, i ); } inline DualNumber operator/(const DualNumber &w1, const DualNumber &w2) // w/q { double r=w1.real()/w2.real(); double i=(-w1.real()*w2.infi()+w1.infi()*w2.real() )/( w2.real()*w2.real() ) ; return DualNumber( r, i ); } inline DualNumber operator/(const DualNumber &w, const double &a) // w/a { double r=w.real()/a; double i=w.infi()/a; return DualNumber( r, i ); } inline DualNumber operator/(const double &a,const DualNumber &w) // a/w { double r=a/w.real(); double i=-a*w.infi()/(w.real()*w.real()); return DualNumber( r, i ); } /* We also need to define the meaning of mathematical functions that take dual numbers and return dual number. Remember, if w=x+ey then f(w) = f(x) + f'(x)ye by Taylor expansion Below I define exp(w) and sqrt(w). This means we can do everything with +,-,*,/ exp(w) and sqrt(w) if you want to use other functions such as e.g. sin(w), cos(w), log(w), pow(w,a), you need to implement their Dual number meaning. Otherwise yoǘ'll get errors like error: no matching function for call to ‘sin(DualNumber&)’ */ inline DualNumber exp(const DualNumber &w) { double r=exp( w.real()) ; double i=exp( w.real())*w.infi() ; // dexp(x)/dx = exp(x) return DualNumber( r, i ); } inline DualNumber sqrt(const DualNumber &w) { double r=sqrt( w.real()) ; double i=0.5/sqrt( w.real()) *w.infi() ; // d x^(1/2)/dx = 0.5/sqrt(x) return DualNumber( r, i ); } /* DualNumber sin(const DualNumber &w) { double r= ..... ; double i= ..... ; return DualNumber( r, i ); } : /* /* Finally to output dual numbers as other numbers, we need to define the stream output operator << for dual numbers, which also means we can write these. */ inline ostream& operator<<(ostream &o, const DualNumber &w) { o << w.real() << "+" << w.infi() << "e"; return o; }