How to solve a system of nonlinear factorial equations in Python

Question

equations

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

Solve for the unknowns n and c!

Answer 1

At a first look, some information is missing regarding how your application defines an acceptable solution, but I am going to make some assumptions and provide a (somewhat naive) approximate answer.

Context first: the problem, posed as is, is a set of 2 algebraic equations in 2 unknowns, to be solved on the 2D grid of non-negative integers. I am assuming you meant the unknowns are non-negative integers, and not, say, real numbers, from the problem formulation, involving the combinatorial terms. Further, note that whenever n is smaller than c, the summation will eventually involve factorial operations on a negative number, formally undefined, hence I will further suggest assuming that we solve over the triangular domain of non-negative integer pairs where n is greater than or equal to c.

As such, to refine the definition of solution, one would have to specify an allowed numeric tolerance (requiring exact equality or allowing, from your application point of view, some error). Unless there is some guarantee, from the problem domain, of the existence of integers solving exactly, I am proceeding to assume that an approximation with known error is meaningful/valuable. Finally: I am going to assume you are looking for A solution (not all solutions) and I am going to limit the search space, as you will see in the simple script below.

With the above said, I don't have any algebraic manipulation to suggest that would lead to a solution, and am not aware of a commercial/open source solver that can solve this non-linear problem over an integer grid. Surprisingly, choosing a reasonable grid and just exhaustively checking some options, yields the following approximate solution:

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

This of course does not claim that there cannot be more accurate approximations for some part of the grid not visited. Such guarantees, if at all, can be achieved with further analysis, such as monotone behavior if the right hand side of the equations w.r.t n and c, etc.

To arrive at the simple approximation suggested above, I treated your right hand sides as functions of n and c, and performed some exhaustive evaluations, keeping track of the error w.r.t to the left hand side required value. The error is then defined by the L2 norm of the error vector - simply the difference between the RHS evaluated at the current location, and the LHS of each of the equations. This is done by the following script, which of course is not scalable for large values of n and c, and the fact that the closest candidate was in the interior of the search space, for all I know, was just chance:

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, c, p1)
    eq2_val = eval_eq(n, c, 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}')