C#凹凸曲线求拐点算法

  1         public static double diff(Func<double, double> f, double x, double h)
  2         {
  3             return (f(x+h)-f(x-h))/(2*h);
  4         }
  5 
  6         public static double diff2(Func<double, double> f, double x, double h)
  7         {
  8             return (f(x+h)-2*f(x)+f(x-h))/(h*h);
  9         }
 10 
 11 
 12         enum PointType { None, Minimum, Maximum, Inflection };
 13 
 14         private static PointType getPointType(Func<double, double> f, double x, double h)
 15         {
 16             double f1 = diff(f, x, h);
 17             double f2 = diff2(f, x, h);
 18 
 19             if (Math.Abs(f2) < 1e-6)
 20             {
 21                 return PointType.None;
 22             }
 23             else if (f2 > 0)
 24             {
 25                 return PointType.Minimum;
 26             }
 27             else if (f2 < 0)
 28             {
 29                 return PointType.Maximum;
 30             }
 31             else if (f1 != 0)
 32             {
 33                 return PointType.Inflection;
 34             }
 35             else
 36             {
 37                 return PointType.None;
 38             }
 39         }
 40 
 41         // List<double> findTurningPoints = null;
 42         //Func<double, double> t = x => x * x * x - 3 * x * x + 2 * x + 1;
 43 
 44         //List<double> points = findTurningPoints(x => (x * x * x - 3 * x * x + 2 * x + 1),-10,10,0.01);
 45 
 46 
 47 
 48 
 49         public static List<double> findTurningPoints(Func<double, double> f, double a, double b, double h)
 50         {
 51             List<double> points = new List<double>();
 52 
 53             double fa = f(a);
 54             double fb = f(b);
 55 
 56             if (fa == fb)
 57             {
 58                 return points;
 59             }
 60 
 61             if (fa < fb)
 62             {
 63                 for (double x = a; x < b; x += h)
 64                 {
 65                     double fx = f(x);
 66                     double fxh = f(x + h);
 67 
 68                     if (fx < fxh)
 69                     {
 70                         continue;
 71                     }
 72 
 73                     double x1 = x - h;
 74                     double x2 = x + h;
 75 
 76                     while (x2 - x1 > 1e-6)
 77                     {
 78                         double xm = (x1 + x2) / 2;
 79                         double f1 = f(xm - h);
 80                         double f2 = f(xm);
 81                         double f3 = f(xm + h);
 82 
 83                         if (f2 < f1 && f2 < f3)
 84                         {
 85                             points.Add(xm);
 86                             break;
 87                         }
 88                         else if (f1 < f3)
 89                         {
 90                             x2 = xm;
 91                         }
 92                         else
 93                         {
 94                             x1 = xm;
 95                         }
 96                     }
 97                 }
 98             }
 99             else
100             {
101                 for (double x = a; x < b; x += h)
102                 {
103                     double fx = f(x);
104                     double fxh = f(x + h);
105 
106                     if (fx > fxh)
107                     {
108                         continue;
109                     }
110 
111                     double x1 = x - h;
112                     double x2 = x + h;
113 
114                     while (x2 - x1 > 1e-6)
115                     {
116                         double xm = (x1 + x2) / 2;
117                         double f1 = f(xm - h);
118                         double f2 = f(xm);
119                         double f3 = f(xm + h);
120 
121                         if (f2 > f1 && f2 > f3)
122                         {
123                             points.Add(xm);
124                             break;
125                         }
126                         else if (f1 > f3)
127                         {
128                             x2 = xm;
129                         }
130                         else
131                         {
132                             x1 = xm;
133                         }
134                     }
135                 }
136             }
137 
138             for (int i = 0; i < points.Count; i++)
139             {
140                 PointType type = getPointType(f, points[i], h);
141 
142                 if (type == PointType.None)
143                 {
144                     points.RemoveAt(i);
145                     i--;
146                 }
147             }
148 
149             return points;
150         }