# Approximating solutions to quadratic recurrence boundary value problem

I have a branching process problem that has been reduced to solving a system of recurrences defined for $0 \leq i \leq T$, for some large positive integer $T$ and fixed constants $d, \mu, \beta$ by

$$p_{i} = d + (1-d)((1-\mu-\beta)p_{i}^{2} + \mu p_{i}p_{i+1} + \beta p_{i}p_{i-1}), p_{-1} = 0, p_{T} = 1$$

Here, $0 < d < 1$ and $0 < \mu + \beta < 1$. Since this is a non-homogeneous quadratic recurrence, I highly doubt that there exist closed form solutions for the $p_i$. Does anyone know of a way to instead approximate the $p_i$, either with an analytical approach or with a computer algorithm? I see that there is a post here regarding a more general form of this problem: Algorithm for solving system of quadratic equations and linear equations. However, given the extreme sparsity of the corresponding matrix here, I feel like there must be an easier solution. Any help is appreciated. Thanks.

• What have you tried so far? If you just try a Newton-type method, from any nonlinear optimization library, does it fail? When I tried solving it like that in Mathematica, it worked and it seems mostly constant $\frac{d}{1-d}$, with two boundary layers at $0,T$, so perturbation theory might work too. Sep 9 '17 at 18:28
• I edited the title to make it clearer that it's a BVP, which is where the difficulties might come from. I think more often than not, "solving a recurrence" means solving an IVP. Sep 9 '17 at 18:46
• Thanks for the response. I haven't really tried anything significant so far - I know very little in the ways of solving polynomial systems. I was thinking of using Grobner bases, but these equations seem like they would hit the exponential worst case for that method. Can you go a bit further into detail about the Newton-type methods you're describing?
– Alex
Sep 9 '17 at 19:06

Newton-type methods are standard methods for solving system of nonlinear equations numerically. To make sure they can handle this, I tried solving your problem with NLsolve.jl in julia, and it seems to work just fine:

using NLsolve
using LineSearches
using ForwardDiff

function solve(d::Real, μ::Real, β::Real, T::Integer)
@assert 0 < d < 1
@assert 0 < μ + β < 1
n = T + 1
rhs(p₋, p, p₊) = -p + d + (1-d)*p*((1-μ-β)*p+μ*p₊+β*p₋)
function f!(p, val)
val = p
for k = 2:T
val[k] = rhs(p[k-1], p[k], p[k+1])
end
val[n] = p[n] - 1
end
function g!(p, J)
J[1,1] = J[n,n] = 1
for k = 2:T
J[k,k-1:k+1] = ForwardDiff.gradient(x -> rhs(x...), p[k-1:k+1])
end
end
Jₚ = spzeros(n, n)
Jₚ[1,1] = Jₚ[n,n] = 1
for k = 2:T; Jₚ[k,k-1:k+1] = 1; end
p₀ = ones(n) * min(1, d/(1-d))
df = DifferentiableGivenSparseMultivariateFunction(f!, g!, Jₚ)
@time sol = nlsolve(df, p₀, method=:newton, linesearch! = LineSearches.strongwolfe!, show_trace = true)
g!(sol.zero, Jₚ)
@printf "Jacobian condition number: %.2e\n" cond(Array(Jₚ), Inf)

sol
end


I noticed that the case $\mu+\beta\approx 1$ seems to be the hard case here, and the solution depends on a good initial guess, because the equation is stiff and there are boundary layers at $0$ and $T$. So this isn't a complete answer, it just works in the non-extreme cases.

Generally speaking, one can't expect such a problem to have a closed form solution. Also, Groebner bases are feasible on relatively small systems, not like this one. It is true that the Jacobian here is sparse, but that's not really enough. Simply having a sparse system of nonlinear equations doesn't necessarily mean there will be a simple efficient way of solving them.

Also, as mentioned in this question on MO, multivariate systems of quadratic equations are a very difficult type of problem in general, at least NP-hard and possibly harder.

• Thanks so much! Just wondering, would the equivalent of the nlsolve Julia method in Python be the fsolve method, if you are familiar with Python?
– Alex
Sep 10 '17 at 20:12
• @Alex My advice would be to just try directly all the methods at docs.scipy.org/doc/scipy/reference/optimize.html, it is easy enough, and see which ones work well enough. I'm not sure what the direct equivalent would be, but you don't need a direct equivalent anyway. Sep 10 '17 at 20:30

In addition to Kirill's excellent answer, one thing I think work mentioning is Smale's $\alpha$-theory. One of the hard parts about applying Newton's method is that, if you don't make a really good initial guess, the iteration can diverge. To remedy this problem, it's common to use a damped Newton iteration, but to select the damping parameter all you've got are heuristics.

Smale's theory tells you sufficient conditions for when an initial guess to Newton's method is in the quadratic convergence region of some root of a nonlinear equation when the function is analytic. The sufficiency condition requires that you have uniform bounds on the Taylor coefficients of the function in question. For polynomials, the conditions are especially easy to evaluate since the Taylor coefficients are 0 after the degree of the polynomial, in your case 2. Additionally, the sparsity of the system you're solving makes evaluating the conditions to apply Smale's theory even simpler.

My old officemate has some slides that I think are a really good reference for the subject, and he also has some example code in python.