auto_differentiate.cpp
#include <math.h> | |
#include <stdio.h> | |
template <class Number> | |
Number f(Number x) { | |
Number z = x; | |
for (int i = 0; i < 4; ++i) { | |
z = z*z - 0.1 + sin(z*3); | |
} | |
if (z < 1) { | |
return z; | |
} else { | |
return z / (z + 1); | |
} | |
} | |
struct RealDiff { | |
double x; | |
double dx; | |
RealDiff(double x, double dx) : x(x), dx(dx) {} | |
typedef RealDiff Me; | |
Me operator-(double y) { return Me(x - y, dx); } | |
Me operator+(double y) { return Me(x + y, dx); } | |
Me operator*(double y) { return Me(x * y, dx * y); } | |
Me operator/(double y) { return Me(x / y, dx / y); } | |
Me operator+(Me y) { return Me(x + y.x, dx + y.dx); } | |
Me operator*(Me y) { return Me(x * y.x, dx*y.x + y.dx*x); } | |
Me operator/(Me y) { return Me(x / y.x, (dx*y.x - y.dx*x) / (y.x*y.x)); } | |
bool operator<(double y) { return x < y; } | |
friend Me operator/(double x, Me y) { return Me(x,0)/y; } | |
friend Me sin(Me y) { return Me(sin(y.x), cos(y.x)*y.dx); } | |
}; | |
int main() { | |
double a = -3, b = +3; | |
for (int i = 0; i < 10000; ++i) { | |
double t = i / (10000.0 - 1); | |
double x = a + (b-a)*t; | |
RealDiff r = f(RealDiff(x, 1.0)); | |
double numeric_diff = (f(x + 1e-5) - f(x)) / 1e-5; | |
printf("%g %g %g %g\n", x, r.x, r.dx, numeric_diff); | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment