// Redesign and adaptation by Babak Sheykhan // Begin : May 12, 2018 2:35 PM // End : June 18, 2018 1:25 PM #define _USE_MATH_DEFINES #include #include #include using namespace std; #define C(N) char(N) /*-----------------------------*/ /* */ /* Background Visual Elements */ /* */ /*-----------------------------*/ #define NAME_1 C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), 'A', 'r', 'c', 'h', 'a', 'i', 'c' #define NAME_2 'P', 'l', 'o', 't', 't', 'e', 'r', C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219) #define X_SCALE(N, F) C(219), 'X', 's', 'c', 'a', 'l', 'e', ' ', N + '0', '.', F + '0', 'x', C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219) #define Y_SCALE(N, F) C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), 'Y', 's', 'c', 'a', 'l', 'e', ' ', N + '0', '.', F + '0', 'x', C(219) #define EDGE C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219), C(219) #define GRID_A_L C(219), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196) #define GRID_B_L C(219), C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ' #define X_AXIS_L C(219), C(216), C(205), C(205), C(205), C(216), C(205), C(205), C(205), C(216), C(205), C(205), C(205), C(216), C(205), C(205), C(205), C(216), C(205), C(205), C(205) #define X_AXIS_R C(205), C(205), C(205), '1', C(205), C(205), C(205), '2', C(205), C(205), C(205), '3', C(205), C(205), C(205), '4', C(205), C(205), C(205), '5', C(219) #define GRID_A_R C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(196), C(196), C(196), C(197), C(219) #define GRID_B_R ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), ' ', ' ', ' ', C(179), C(219) /* .................................... ... ____________________________ ... ...| CurveArea (UML) |... ...+----------------------------+... ...| [-] subIntervals : Integer |... ...| <> |... ...|____________________________|... ...| [+] calculateArea( |... ...| Double f(Double, Boolean),|... ...| Double f(Double, Double), |... ...| Double, |... ...| Double, |... ...| Integer, |... ...| Double, |... ...| Float, |... ...| Float |... ...| ) <> |... ...| |... ...| [-] plotter( |... ...| Double f(Double, Boolean),|... ...| Float, |... ...| Float |... ...| ) <> |... ...| |... ...| [-] sliceArea( |... ...| Double f(Double, Boolean),|... ...| Double, |... ...| Double, |... ...| Double |... ...| ) : Double <> |... ...| |... ...| [-] calcErr(Double, Double)|... ...| : Double <> |... ...|____________________________|... .................................... */ /*-----------------------------*/ /* */ /* CurveArea.h */ /* */ /*-----------------------------*/ class CurveArea { public: // Calculate and compare the area under the curve using 4 different approaches. static void calculateArea( double (*f_est)(double, bool), // Estimate function double (*f_act)(double, double), // Actual function double, // x-axis upper range double, // x-axis lower range int, // Sampling rate (interval, partition, slice, or whatever you call it) double, // Maximum error percentage of the slice float, // x-axis zoom (scale) float // y-axis zoom (scale) ); private: static void plotter( double (*f)(double, bool), // Estimate function float, // x-axis zoom (scale) float // y-axis zoom (scale) ); static double sliceArea( double (*f)(double, bool), // Estimate function double, // Slice's x1 value double, // Slice's x2 value double // Maximum error percentage of the slice ); static double calcErr(double, double); static int subIntervals; }; /*-----------------------------*/ /* */ /* CurveArea.cpp */ /* */ /*-----------------------------*/ int CurveArea::subIntervals = 0; void CurveArea::calculateArea( double (*f_est)(double, bool), double (*f_act)(double, double), double xMin, double xMax, int interval, double maxSliceErr, float xZoom, float yZoom ) { plotter(f_est, xZoom, yZoom); // calculate the width of a single rectangle (delta x). double dx = (xMax - xMin) / interval; // Add up the triangles', trapezoids', and adaptive adaptive mid-point areas. double total_rectangular_area = 0.0; double total_trapezoidal_area = 0.0; double total_adaptiveMid_area = 0.0; double total_actual_area = f_act(xMin, xMax); double x = xMin; maxSliceErr /= 100.0; for (int i = 1; i <= interval; ++i) { // Riemann Sum // https://en.wikipedia.org/wiki/Riemann_sum // // Using Left Riemann sum using the rectangle and trapezoid rules. // How does it work? Well, it uses a series of rectangles, trapezoids etc. of uniform width to approximate the area under a curve. // After calculating the height of the rectangle by calling f_est and multiplying it by dx (width), it's going to add up to total area. total_rectangular_area += dx * f_est(x, false); // Trapezoid area // // .-----a----. // /| \ // / h \ // /__|_____b_______\ // // a + b // area = ---------- * h // 2 // Transforming the above formula to what we are after gives us: // h = dx // a = f(x) // b = f(x + 1) // // f(x) + f(x + 1) // area = ----------------- * dx // 2 total_trapezoidal_area += (f_est(x, false) + f_est(x + dx, false)) / 2.0 * dx; total_adaptiveMid_area += sliceArea(f_est, x, x + dx, maxSliceErr); x += dx; } double rectangle_err = calcErr(total_rectangular_area, total_actual_area); double trapezoid_err = calcErr(total_trapezoidal_area, total_actual_area); double aMidpoint_err = calcErr(total_adaptiveMid_area, total_actual_area); cout << "Input range [" << xMin << ", " << xMax << "]\n" << setprecision(4); cout << "+-----------+---------+--------+-------+-----+\n"; cout << "| method |initSlice|finSlice| area |Err %|\n"; cout << "+-----------+---------+--------+-------+-----+\n"; cout << "| Rectangle |" << setw(5) << interval << setw(10) << "--" << setw(10) << total_rectangular_area << setw(7) << rectangle_err << endl; cout << "+-----------+\n"; cout << "| Trapezoid |" << setw(5) << interval << setw(10) << "--" << setw(10) << total_trapezoidal_area << setw(7) << trapezoid_err << endl; cout << "+-----------+\n"; cout << "|AdaptiveMid|" << setw(5) << interval << setw(10) << subIntervals << setw(10) << total_adaptiveMid_area << setw(9) << aMidpoint_err << endl; cout << "+-----------+\n"; cout << "| Calculus |" << setw(5) << "--" << setw(10) << "--" << setw(10) << total_actual_area << setw(7) << "0.0" << endl; cout << "+-----------+--------------------------------+\n\n"; } void CurveArea::plotter( double (*f)(double, bool), float xScale, float yScale ) { // reset to default values if scales was out of bound xScale = (xScale >= 10.0f || xScale <= 0.0f) ? 1.0f : xScale; yScale = (yScale >= 10.0f || yScale <= 0.0f) ? 1.0f : yScale; int xScale_frac = static_cast((xScale - static_cast(xScale)) * 10.0f); int yScale_frac = static_cast((yScale - static_cast(yScale)) * 10.0f); const unsigned short xAxis = 43; const unsigned short yAxis = 35; char background[yAxis][xAxis] = { {EDGE, C(219), EDGE}, {NAME_1, ' ', NAME_2}, {EDGE, C(219), EDGE}, {X_SCALE(int(xScale), xScale_frac), C(219), Y_SCALE(int(yScale), yScale_frac)}, {EDGE, C(219), EDGE}, {GRID_A_L, '7', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, '6', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, '5', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, '4', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, '3', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, '2', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, '1', GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {X_AXIS_L, C(206), X_AXIS_R}, // X-Axis {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {GRID_B_L, C(186), GRID_B_R}, {GRID_A_L, C(215), GRID_A_R}, {EDGE, C(219), EDGE} }; // Preparing the plot // Starting from -5 to 5 with 0.25 increment float x = -5.0f; for (int i = 1; i < (xAxis - 1); ++i) { double fx = f(x * xScale, true) * yScale; double fraction = fx - static_cast(fx); // Rounding to the near 0.5 // // Is the fraction between [0.25 and 0.75]? // If yes, then take a floor of it and then add 0.5. Othrwise add 0.5 to it then take a floor. fx = ( abs(fraction) < 0.75f && abs(fraction) > 0.25f ? floor(fx) + 0.5f : floor(fx + 0.5f) ); // Starting from 7 to -7 with 0.5 increment double y = 7.0f; for (int j = 5; j < (yAxis - 1); ++j) { if (y == fx) background[j][i] = '*'; // walking pace on y-axis (up to down) y -= 0.5; } // Walking pace on the x-axis (left to right) x += 0.25; } // plot -> console for (int i = 0; i < yAxis; ++i) { for (int j = 0; j < xAxis; ++j) cout << background[i][j]; cout << endl; } cout << endl; } double CurveArea::sliceArea( double (*f)(double, bool), double x1, double x2, double maxSliceErr ) { // Calculate the function at the endpoints and the midpoint. double y1 = f(x1, false); double y2 = f(x2, false); double xMid = (x2 + x1) / 2.0; double yMid = f(xMid, false); // Calculate the area for the large slice and two subslices using trapezoid formula. double area12 = (y1 + y2) / 2.0 * (x2 - x1); double area1mid = (y1 + yMid) / 2.0 * (xMid - x1); double areaMid2 = (yMid + y2) / 2.0 * (x2 - xMid); double area1mid2 = area1mid + areaMid2; // See how close we are to actual value double err = (area1mid2 - area12) / area12; // See if error percentage is small enough if (abs(err) < maxSliceErr) { subIntervals += 2; return area1mid2; } return sliceArea(f, x1, xMid, maxSliceErr) + sliceArea(f, xMid, x2, maxSliceErr); } double CurveArea::calcErr(double est_area, double act_area) { return (est_area - act_area) / act_area * 100.0; } /* .................................... ... ____________________________ ... ...| CurveArea (UML) |... ...+----------------------------+... ...| [No Attribute] |... ...|____________________________|... ...| [+] f1(Double, Boolean = F)|... ...| <> |... ...| [+] f2(Double, Boolean = F)|... ...| <> |... ...| [+] f3(Double, Boolean = F)|... ...| <> |... ...| |... ...| [+] actualF1( |... ...| Double, Double) <> |... ...| |... ...| [+] actualF2( |... ...| Double, Double) <> |... ...| |... ...| [+] actualF3( |... ...| Double, Double) <> |... ...|____________________________|... .................................... */ /*-----------------------------*/ /* */ /* F.h */ /* */ /*-----------------------------*/ class F { public: // Candidate functions for estimating the actual areas static double f1(double, bool = false); static double f2(double, bool = false); static double f3(double, bool = false); // Integral of candidate functions for calculating the excact areas static double actualF1(double, double); static double actualF2(double, double); static double actualF3(double, double); }; /*-----------------------------*/ /* */ /* F.cpp */ /* */ /*-----------------------------*/ double F::f1(double in, bool isPlotterCalled) { // need flag double out = sin(in); return isPlotterCalled ? out : abs( out ); } double F::f2(double in, bool isPlotterCalled) { // need flag double out = in * (in - 1.0) * (in - 2.0); // in^3 - 3*in^2 + 2*in return isPlotterCalled ? out : abs( out ); } double F::f3(double in, bool isPlotterCalled) { // don't need flag double out = 1.0 + in + sin(2.0 * in); return out; } double F::actualF1(double xMin, double xMax) { if (xMin < -M_PI || xMax > M_PI) { cout << "Pre-defined range for actual area is [-Pi - Pi]\n"; xMin = -M_PI; xMax = M_PI; } double xMid = 0.0; double positive_peak = abs( (-cos(xMid)) - (-cos(xMin)) ); double negative_peak = abs( (-cos(xMax)) - (-cos(xMid)) ); double out = positive_peak + negative_peak; return out; } double F::actualF2(double xMin, double xMax) { if (xMin < 0.0 || xMax > 2.0) { cout << "Pre-defined range for actual area is [0 - 2]\n"; xMin = 0.0; xMax = 2.0; } double xMid = 1.0; double positive_peak = abs( ((xMid * xMid * xMid * xMid) / 4.0 - (xMid * xMid* xMid) + (xMid * xMid)) ); double negative_peak = abs( ((xMax * xMax * xMax * xMax) / 4.0 - (xMax * xMax* xMax) + (xMax * xMax)) - ((xMid * xMid * xMid * xMid) / 4.0 - (xMid * xMid* xMid) + (xMid * xMid)) ); double out = positive_peak + negative_peak; return out; } double F::actualF3(double xMin, double xMax) { double out = ( xMax + (xMax * xMax) / 2.0 - cos(2.0 * xMax) / 2.0 ) - ( xMin + (xMin * xMin) / 2.0 - cos(2.0 * xMin) / 2.0 ); return out; } /*-----------------------------*/ /* */ /* Client */ /* */ /*-----------------------------*/ int main() { const int sample_rate = 2; const double slice_err = 1.0; CurveArea::calculateArea( F::f3, F::actualF3, 0.0, 5.0, sample_rate, slice_err, 1.1f, 1.2f ); CurveArea::calculateArea( F::f2, F::actualF2, 0.0, 0.75, sample_rate, slice_err, 0.5f, 2.5f ); /*CurveArea::calculateArea( F::f1, F::actualF1, 0, M_PI / 2.0, sample_rate, slice_err, 0.5f, 2.5f );*/ } /* https://en.wikipedia.org/wiki/Riemann_sum https://en.wikipedia.org/wiki/Adaptive_quadrature */

Area under the curve