寻找给定区间内黑盒函数的首个根与局部极值的高效方法

针对你在大区间内寻找组合函数首个根和极值的需求，结合函数是常函数+已知周期的周期函数这一关键特性，我们可以跳出常规根查找的局限，利用周期性大幅压缩搜索范围，完全能达到毫秒级的运行要求。下面是具体的解决方案：

你的函数可以表示为：f(x) = C + P(x)，其中C是常函数，P(x)是已知周期的周期函数集合的组合。这类函数的根、极值点都会随周期重复出现——因为P(x)的周期性会直接传递给f(x)和它的导数f’(x)。

所以我们不需要在整个[0, 1e8]区间搜索，只需要：

1. 计算所有已知周期的最小公倍数（LCM），得到函数的基本周期L；
2. 在单个周期[0, L]内找出所有根和极值点；
3. 把周期内的结果映射到大区间，直接取首个出现的点即可。

这个思路能把搜索范围从1e8压缩到L（通常L远小于1e8），效率提升几个数量级。

步骤1：在单个周期内枚举所有根 #

因为周期L很小，我们可以用Brent's方法（比二分法快得多，结合了割线法和反二次插值，适合单根/多根场景）在[0, L]内找出所有根。具体操作：

• 先把[0, L]分割成若干小片段，检测每个片段内的函数值符号变化（符号变化说明存在根）；
• 对有符号变化的片段用Brent's方法精确定位根。

步骤2：映射到大区间找首个根 #

周期内的根会以x_i + n*L（n为非负整数）的形式在大区间重复出现。因此：

• 如果大区间从0开始，首个根就是周期内最小的那个根；
• 如果大区间从start > 0开始，计算start所在的周期位置，在[start, start+L]内找最小的根即可。

特殊情况：无理数周期 #

如果周期是无理数（无法计算精确LCM），可以先在[0, K*T]（K取10~20，T为其中一个周期）内搜索，再验证后续区间的重复性；或者用数值方法计算近似周期后再缩小范围。

极值点对应导数f’(x)=0的点，而f’(x)同样是周期函数，方法和找根类似：

1. 计算数值导数：因为是黑盒函数，用中心差分近似导数：f’(x) = (f(x+h) - f(x-h))/(2h)，h取1.0即可（符合你的精度要求）；
2. 在单个周期内找导数的根：用Brent's方法在[0, L]内找出所有导数为0的点；
3. 判断极值类型：通过比较候选点邻域的函数值（或二阶数值导数），区分极小/极大值；
4. 映射到大区间：取周期内最小的极小/极大值点，就是大区间内的首个极值点。

为了确保算法稳定高效，需要补充几个前提：

• 周期的准确性：必须确认所有已知周期确实是P(x)的周期，可以通过验证|f(x) - f(x+T)| < 1e-6（远低于你的精度要求）来确认；
• 函数的连续性与可微性：你已经提到这一点，这是数值方法的基础；
• 区间起始点明确：如果大区间不是从0开始，需要明确起始位置，调整搜索的周期区间；
• 精度阈值设置：你的精度要求是1.0，所以数值计算时把收敛阈值设为0.5即可，能大幅加快收敛速度。

下面是核心逻辑的Java实现，你可以根据自己的黑盒函数调整：

import java.util.ArrayList;
import java.util.List;
import java.util.function.Function;
import java.util.stream.Collectors;

public class PeriodicFunctionSolver {
    // 计算整数周期的LCM
    private static long lcm(long a, long b) {
        return a * b / gcd(a, b);
    }

    private static long gcd(long a, long b) {
        while (b != 0) {
            long temp = b;
            b = a % b;
            a = temp;
        }
        return a;
    }

    // 计算多个整数周期的基本周期
    public static double computeBasicPeriod(List<Long> periods) {
        long lcm = 1;
        for (long p : periods) {
            lcm = lcm(lcm, p);
        }
        return (double) lcm;
    }

    // Brent's方法找区间内所有根
    public static List<Double> findAllRoots(double a, double b, Function<Double, Double> f, double tolerance) {
        List<Double> roots = new ArrayList<>();
        double step = (b - a) / 500; // 根据周期大小调整步长
        double x = a;
        double fx = f.apply(x);

        while (x <= b) {
            double nextX = Math.min(x + step, b);
            double fNextX = f.apply(nextX);
            if (fx * fNextX <= 0) {
                double root = brentSingleRoot(x, nextX, f, tolerance);
                if (!Double.isNaN(root)) roots.add(root);
            }
            x = nextX;
            fx = fNextX;
        }
        // 去重（避免相邻片段找到同一个根）
        return roots.stream().distinct().collect(Collectors.toList());
    }

    // Brent's方法找单个根
    private static double brentSingleRoot(double a, double b, Function<Double, Double> f, double tolerance) {
        double fa = f.apply(a);
        double fb = f.apply(b);
        if (fa * fb > 0) return Double.NaN;

        if (Math.abs(fa) < Math.abs(fb)) {
            double temp = a; a = b; b = temp;
            double tempF = fa; fa = fb; fb = tempF;
        }

        double c = a, fc = fa;
        boolean mflag = true;
        double d = 0;

        while (Math.abs(b - a) > tolerance && fb != 0) {
            double x;
            if (fa != fc && fb != fc) {
                // 反二次插值
                x = a*fb*fc/((fa-fb)*(fa-fc)) + b*fa*fc/((fb-fa)*(fb-fc)) + c*fa*fb/((fc-fa)*(fc-fb));
            } else {
                // 割线法
                x = b - fb*(b-a)/(fb-fa);
            }

            // 确保x在合理范围内，否则退化为二分法
            if (x < Math.min(a,b) || x > Math.max(a,b) || 
                (mflag && Math.abs(x-b) >= Math.abs(b-c)/2) || 
                (!mflag && Math.abs(x-b) >= Math.abs(c-d)/2)) {
                x = (a + b)/2;
                mflag = true;
            } else {
                mflag = false;
            }

            double fx = f.apply(x);
            d = c; c = b; fc = fb;
            if (fa * fx < 0) { b = x; fb = fx; }
            else { a = x; fa = fx; }

            if (Math.abs(fa) < Math.abs(fb)) {
                double temp = a; a = b; b = temp;
                double tempF = fa; fa = fb; fb = tempF;
            }
        }
        return b;
    }

    // 计算数值导数
    public static Function<Double, Double> getDerivative(Function<Double, Double> f, double h) {
        return x -> (f.apply(x+h) - f.apply(x-h))/(2*h);
    }

    // 判断是否为局部极小值
    public static boolean isLocalMin(double x, Function<Double, Double> f, double h) {
        double fx = f.apply(x);
        return f.apply(x-h) > fx && f.apply(x+h) > fx;
    }

    // 判断是否为局部极大值
    public static boolean isLocalMax(double x, Function<Double, Double> f, double h) {
        double fx = f.apply(x);
        return f.apply(x-h) < fx && f.apply(x+h) < fx;
    }

    public static void main(String[] args) {
        // 示例：f(x) = 3 + sin(x) + cos(2x)，周期为2π≈6.283
        List<Long> periods = List.of(2L, 1L); // 模拟sin(2πx)和cos(4πx)的周期
        double basicPeriod = computeBasicPeriod(periods);
        Function<Double, Double> f = x -> 3 + Math.sin(x) + Math.cos(2*x);

        // 找首个根
        List<Double> roots = findAllRoots(0, basicPeriod, f, 0.5);
        roots.stream().min(Double::compare)
            .ifPresent(root -> System.out.println("首个根：" + Math.round(root)));

        // 找首个局部极小值
        Function<Double, Double> df = getDerivative(f, 1.0);
        List<Double> criticalPoints = findAllRoots(0, basicPeriod, df, 0.5);
        criticalPoints.stream()
            .filter(x -> isLocalMin(x, f, 1.0))
            .min(Double::compare)
            .ifPresent(min -> System.out.println("首个局部极小值点：" + Math.round(min)));
    }
}

内容的提问来源于stack exchange，提问作者Abastro