如何在 Python 中求解非线性阶乘方程组

问题描述

alpha=0.05,beta=0.1,p1=0.015,p2=0.1

求解未知数 n 和 c！

解决方法

乍一看，缺少关于您的应用程序如何定义可接受解决方案的一些信息，但我将做出一些假设并提供（有点幼稚的）近似答案。

上下文优先：问题按原样提出，是一组 2 个未知数中的 2 个代数方程，要在非负整数的 2D 网格上求解。我假设您的意思是未知数是非负整数，而不是来自问题公式的实数，涉及组合项。此外，请注意，当 n 小于 c 时，求和最终将涉及对正式未定义的负数进行阶乘运算，因此我将进一步建议假设我们在非负整数对的三角域上求解，其中 n 大于或等于 c。

因此，要完善解决方案的定义，必须指定允许的数值容差（要求完全相等或从您的应用程序的角度来看允许存在一些错误）。除非从问题域中保证存在精确求解的整数，否则我将继续假设具有已知误差的近似值是有意义的/有价值的。最后：我将假设您正在寻找解决方案（并非所有解决方案）并且我将限制搜索空间，正如您将在下面的简单脚本中看到的那样。

综上所述，我没有任何代数运算来暗示会导致解决方案，并且不知道可以解决整数网格上的这个非线性问题的商业/开源求解器。令人惊讶的是，选择一个合理的网格并详尽地检查一些选项，会产生以下近似解决方案：

Opt n: 52,opt c: 2,min_err (residual L2 norm): 0.0074985711238005955

这当然并不是说对于未访问的网格的某些部分没有更准确的近似值。这种保证，如果有的话，可以通过进一步的分析来实现，例如如果方程的右手边 w.r.t n 和 c 等，则可以实现单调行为。

为了得到上面建议的简单近似值，我将您的右手边视为 n 和 c 的函数，并进行了一些详尽的评估，跟踪误差 w.r.t 到左手边所需的值。然后，误差由误差向量的 L2 范数定义 - 只是在当前位置评估的 RHS 与每个方程的 LHS 之间的差异。这是通过以下脚本完成的，当然，对于 n 和 c 的大值，该脚本不可扩展，而且据我所知，最接近的候选位于搜索空间内部的事实只是偶然：

import math
import numpy as np


def eval_eq(n,c,p):
    res = 0
    for d in range(c + 1):
        res += (math.factorial(n) / (math.factorial(d) * (math.factorial(n - d)))) * (p ** d) * ((1. - p) ** (n - d))
    return res


def eval_err_norm(n,c):
    p1 = 0.015
    p2 = 0.1
    beta = 0.1
    alpha = 0.05

    eq1_val = eval_eq(n,p1)
    eq2_val = eval_eq(n,p2)
    err1 = eq1_val - (1. - alpha)
    err2 = eq2_val - beta
    l2_err = np.linalg.norm(np.array([err1,err2]))
    return l2_err


if __name__ == '__main__':
    err = np.inf
    opt_n = None
    opt_c = None

    for n_ in range(100):
        for c_ in range(100):
            if n_ < c_:
                continue
            else:
                cur_err = eval_err_norm(n_,c_)
                if cur_err < err:
                    err = cur_err
                    opt_n = n_
                    opt_c = c_
    print(f'Opt n: {opt_n},opt c: {opt_c},min_err (residual L2 norm): {err}')

factorial nonlinear-equation python