java - 实现自适应函数绘图

Question

我正在尝试使用这两个示例中的伪代码来实现自适应函数绘图算法(两个示例实际上相同)

https://www.andr.mu/logs/acquiring-samples-to-plot-a-math-function-adaptive/ http://yacas.readthedocs.io/en/latest/book_of_algorithms/basic.html

我遇到的问题(来 self 可能不正确的实现)是:

(1) 正在创建重复的坐标。我知道这是因为 split ，每次 split 的两端都保留，因此一个间隔的结束是另一个间隔的开始 - 相同的 x 值评估两次。但是有没有办法设置算法来避免重复的坐标。我可以避免添加开始或结束坐标作为一个肮脏的修复(请参阅下面代码中的注释)，但随后我会在整个间隔中丢失一个坐标。

(2) 绘图的某些部分缺少本质上对称函数的坐标，这很奇怪。该算法对双方来说应该以相同的方式工作，但事实并非如此。当深度 >= 6 时，这种情况就会开始发生，它对函数的右侧进行额外的分割，而对原点周围的左侧不执行任何操作。远离原点的其余坐标似乎是匹配的。总体而言，右侧似乎比左侧有更多的裂痕。

问题(2)

我的算法实现

static List<Double[]> linePoints;

public static void main(String[] args) {
    linePoints =  new ArrayList<>();

    double startX = -50;
    double endX = 50;

    sampling(startX, endX, depth, tolerance);

    /* Print all points to be plotted - x,y coordinates */
    for (Double[] point : linePoints) {
        System.out.println(point[0]+","+point[1]);
    }
}

/* math function */
public static double f(double x){
    return x*Math.sin(x);
}


static int depth = 6; /* 8 */
static double tolerance = 0.005; /* just a guess */

/* Adaptive sampling algorithm */
/* mostly followed along 2st website and used 1st website variable names  */
public static void sampling(double xa, double xc, int depth, double tolerance){
    /* step (1) of 2nd website - determine mid-intervals */
    double xb = (xa+xc)/2;   /* (xc-xa)/2; tried these from 1st website - didn't work out */
    double xab = (xa+xb)/2;  /* (xb-xa)/2; */
    double xbc = (xb+xc)/2;  /* (xc-xb)/2; */

    /* evaluate the above points using math function - store in array */
    double[] points = new double[5];
    points[0] = f(xa); points[1] = f(xab); points[2] = f(xb); points[3] = f(xbc); points[4] = f(xc);

    /* step (2) of 2nd website */
    if (depth <= 0){
        linePoints.add(new Double[]{xa, points[0]});  /* either I comment out this line for dirty fix */
        linePoints.add(new Double[]{xab, points[1]});
        linePoints.add(new Double[]{xb, points[2]});
        linePoints.add(new Double[]{xbc, points[3]});
        linePoints.add(new Double[]{xc, points[4]});  /* or comment out this line */
    } else {
        /* step (3) of 2nd website */
        int counter = 0;

        for (int i = 1; i < points.length-1; i++){
            /* Check if prev, current, next values are infinite or NaN */
            if (    (Double.isInfinite(points[i-1]) || Double.isNaN(points[i-1])) ||
                    (Double.isInfinite(points[i]) || Double.isNaN(points[i])) ||
                    (Double.isInfinite(points[i+1]) || Double.isNaN(points[i+1]))){
                counter++;
                continue;
            }

            /* Determine the fluctuations - if current is < or > both it's left/right neighbours */
            boolean middleLarger = (points[i] > points[i-1]) && (points[i] > points[i+1]);
            boolean middleSmaller = (points[i] < points[i-1]) && (points[i] < points[i+1]);

            if (middleLarger || middleSmaller){
                counter++;
            }
        }

        if (counter <= 2){  /* at most 2 */
            /* Newton-Cotes quadratures - check if smooth enough */
            double f1 = (3d/8d)*points[0]+(19d/24d)*points[1]-(5d/24d)*points[2]+(1d/24d)*points[3];    /* add 'd' to end of number, otherwise get 0 always */
            double f2 = (5d/12d)*points[2]+(2d/3d)*points[3]-(1d/12d)*points[4];

         if (Math.abs(f1-f2) < tolerance * f2){
                linePoints.add(new Double[]{xa, points[0]});
                linePoints.add(new Double[]{xab, points[1]});
                linePoints.add(new Double[]{xb, points[2]});
                linePoints.add(new Double[]{xbc, points[3]});
                linePoints.add(new Double[]{xc, points[4]});
            } else {
                /* not smooth enough - needs more refinement */
                depth--;
                tolerance *= 2;
                sampling(xa, xb, depth, tolerance);
                sampling(xb, xc, depth, tolerance);
            }

        } else {
            /* else (count > 2), that means further splittings are needed to produce more accurate samples */
            depth--;
            tolerance *= 2;
            sampling(xa, xb, depth, tolerance);
            sampling(xb, xc, depth, tolerance);
        }
    }
}

修复 - 修改我的代码

查看 Gene 的示例，将容差乘以 0.5 而不是 2 似乎可以解决问题 (2)

基因示例是该算法的更好、更清晰的实现，并处理重复的坐标

Answer 1

我认为你已经忠实地实现了，但是算法被破坏了。不对称正交很容易给出不对称的结果。我也有同样的奇怪感觉。

但是您可以通过将重复点维护在排序集中并使用递归分析器运行时端点已被插入的不变量来消除重复点。

下面是更充分地使用现代 Java 功能的简化:

import static java.lang.Double.compare;
import static java.lang.Double.isFinite;
import static java.lang.Math.PI;

import java.util.List;
import java.util.SortedSet;
import java.util.TreeSet;
import java.util.function.DoubleUnaryOperator;
import static java.util.stream.Collectors.toList;
import java.util.stream.DoubleStream;

public class AdaptivePlot {
  private final DoubleUnaryOperator f;
  private final double a;
  private final double c;
  private final SortedSet<Point> plot = new TreeSet<>((s, t) -> compare(s.x, t.x));

  public AdaptivePlot(DoubleUnaryOperator f, double a, double c) {
    this.f = f;
    this.a = a;
    this.c = c;
  }

  public static class Point {
    final double x, y;
    public Point(double x, double y) {
      this.x = x;
      this.y = y;
    }
  }

  public AdaptivePlot computePlot(int depth, double eps) {
    plot.clear();
    Point pa = pointAt(a);
    Point pc = pointAt(c);
    plot.add(pa);
    plot.add(pc);
    computePlot(pa, pc, depth, eps);
    return this;
  }

  public List<Point> getPlot() {
    return plot.stream().collect(toList());
  }

  private Point pointAt(double x) {
    return new Point(x, f.applyAsDouble(x));
  }

  private void computePlot(Point pa, Point pc, int depth, double eps) {
    Point pb = pointAt(0.5 * (pa.x + pc.x));
    Point pa1 = pointAt(0.5 * (pa.x + pb.x));
    Point pb1 = pointAt(0.5 * (pb.x + pc.x));
    plot.add(pb);
    if (depth > 0 && 
        (oscillates(pa.y, pa1.y, pb.y, pb1.y, pc.y) 
            || unsmooth(pa.y, pa1.y, pb.y, pb1.y, pc.y, eps))) {
      computePlot(pa, pb, depth - 1, 2 * eps);
      computePlot(pb, pc, depth - 1, 2 * eps);
    }
    plot.add(pa1);
    plot.add(pb1);
  }

  private static boolean oscillates(
      double ya, double ya1, double yb, double yb1, double yc) {
    return isOscillation(ya, ya1, yb) 
        && isOscillation(ya1, yb, yb1) 
        && isOscillation(yb, yb1, yc);
  }

  private static boolean isOscillation(double ya, double yb, double yc) {
    return !isFinite(ya) || !isFinite(yb) || !isFinite(yc) 
        || (yb > ya && yb > yc) || (yb < ya && yb < yc);
  }

  private static boolean unsmooth(
      double ya, double ya1, double yb, double yb1,double yc, double eps) {
    double y0 = DoubleStream.of(ya, ya1, yb, yb1, yc).min().getAsDouble();
    double [] yg = DoubleStream.of(ya, ya1, yb, yb1, yc).map(y -> y - y0).toArray();
    double q4 = quadrature(yg[0], yg[1], yg[2], yg[3]);
    double q3 = quadrature(yg[2], yg[3], yg[4]);
    return Math.abs(q4 - q3) > eps * q3;
  }

  private static double quadrature(double y0, double y1, double y2, double y3) {
    return 3d/8d * y0 + 19d/24d * y1 - 5d/24d * y2 + 1d/24d * y3;
  }

  private static double quadrature(double y0, double y1, double y2) {
    return 5d/12d * y0 + 2d/3d * y1 - 1d/12d * y2;
  }

  public static void main(String [] args) {
    List<Point> plot = new AdaptivePlot(x -> x * Math.sin(x), -2d * PI, 2d * PI)
        .computePlot(6, 0.005).getPlot();
    for (Point p : plot) {
      System.out.println(p.x + "\t" + p.y);
    }
  }
}