Solve for equation that just above 10000 degree

Question

I am trying to solve a equation whose degree rises as n rises：

n = 10002;
equ = -b;
Do[equ = equ + (1 - b) v^i, {i, n - 1}];
vs = Solve[{equ /. {b -> (n - 2)/n}} == 0, v, Reals];
Do[If[0 <= vs[[j, 1, 2]] <= 1, fans = vs[[j, 1, 2]]], {j, Length[vs]}];
fans = vs

(1)I am looking for the real answers between [0,1];

(2)I don't care about whether is a NSolve or Solve;

(3)I am pretty sure when b=(n-2)/n,there is only 1 answer between [0,1] and it should be very close to 1;

What makes me confused is, When I make n=10001, I can get the answer in 1min, but when I make n=10002, I can't get the answer in 15mins. In my opinion, shouldn't the time be continuous?

What's more, I hope you can give me some advice about improving the efficiency on solving high degree equations, it means a lot to me!

Answer 1

Define the equation exactly, without summing explicitly:

Clear[n, equ];
equ[n_, b_, v_] = -b + (1 - b) Sum[v^i, {i, n - 1}] // FullSimplify
(*    (b - v + v^n - b v^n)/(-1 + v)    *)

Insert the desired value of $b$ and simplify:

equ[n, (n - 2)/n, v] // FullSimplify
(*    (2 + n (-1 + v) - 2 v^n)/(n - n v)    *)

Set the numerator to zero and solve:

vv[n_Integer] := SolveValues[2 + n (-1 + v) - 2 v^n == 0 && 0 < v < 1, v]

This is fast (~millisecond) even for very large values of $n$:

vv[10^6]
(*    Root around 0.999998...    *)

Convert the Root object to a numerical value:

vv[10^6] // N
(*    {0.999998}    *)

There is no closed-form solution that Mathematica can find; but we can find an asymptotic expansion for $n\to\infty$:

w[n_] = 1 + k1 n^-1 + k2 n^-2 + k3 n^-3 /.
  {k1 -> -2 - q,
   k2 -> (q (2 + q)^2)/(2 (1 + q)),
   k3 -> -((q (2 + q)^3 (-2 + q (-1 + 4 q)))/(24 (1 + q)^3))} /. 
  q -> ProductLog[-2/E^2];

See that this solution satisfies the equation asymptotically: (in fact, that's how I found its coefficients $k_j$):

Series[equ[n, (n - 2)/n, w[n]], {n, ∞, 2}] // FullSimplify
(*    O[1/n]^3    *)

Check the values:

vv[100] // N
(*    {0.983977}    *)
w[100] // N
(*    0.983977    *)

Numerical coefficients of the asymptotic expansion:

w[n] // N
(*    1 - 1.59362/n - 0.869277/n^2 - 0.305667/n^3    *)

Update: how to find the asymptotic expansion for $n\to\infty$

Usually, Mathematica offers tools to find asymptotic solutions to difficult equations. Unfortunately, for the given case I cannot make AsymptoticSolve give an asymptotic solution for $v(n)$. SolveAlways wasn't any help either.

But we can do it manually. First, set up the desired asymptotic form:

w[n_] = 1 + k1 n^-1 + k2 n^-2 + k3 n^-3;

Inserting this solution into the equation to be solved (to be set to zero) and series-expanding around $n\to\infty$ gives

Series[equ[n, (n - 2)/n, w[n]], {n, ∞, 2}] // FullSimplify
(*    -((2 - 2 E^k1 + k1)/k1) +
      (2 k2 - E^k1 (k1^3 + 2 k2 - 2 k1 k2))/(k1^2 n) +
      (24 (-k2^2 + k1 k3) + E^k1 (8 k1^5 + 3 k1^6 - 12 k1^3 k2 - 12 k1^4 k2 + 24 k2^2 - 24 k1 (k2^2 + k3) + 12 k1^2 (k2^2 + 2 k3)))/(12 k1^3 n^2) +
      O[1/n]3    *)

We want this expression to be zero for all values of $n$. Setting the first term to zero gives the value $k_1$,

Solve[(-2 + 2 E^k1 - k1)/k1 == 0, k1] // FullSimplify
(*    {{k1 -> -2 - ProductLog[-(2/E^2)]}}    *)

We are lucky here that the ProductLog gave the correct solution branch (see warning).

Insert this value of $k_1$ and repeat the process, to find $k_2$:

w[n_] = 1 + k1 n^-1 + k2 n^-2 + k3 n^-3 /.
        {k1 -> -2 - ProductLog[-(2/E^2)]};
Series[equ[n, (n - 2)/n, w[n]], {n, ∞, 2}] // FullSimplify
(*    ((-2 + E^(2 + ProductLog[-(2/E^2)])) k2 + (2 + ProductLog[-(2/E^2)])^2)/(2 (-1 + E^(2 + ProductLog[-(2/E^2)])) n) + [...]    *)
Solve[((-2 + E^(2 + ProductLog[-(2/E^2)])) k2 + (2 + ProductLog[-(2/E^2)])^2)/(2 (-1 + E^(2 + ProductLog[-(2/E^2)])) n) == 0,
      k2] // FullSimplify
(*    {{k2 -> -((2 + ProductLog[-(2/E^2)])^2/(-2 + E^(2 + ProductLog[-(2/E^2)])))}}    *)

We can continue this process to find $k_3$ and keep going as long as we want to find higher-order coefficients of the asymptotic solution.