Integrate does not converge but NIntegrate does

Question

I have a problem with Mathematica 13.0. The function f(z) does not converge in the range [zcut, 1] analytically but it does work with NIntegrate for fixed values of rho and s, where zcut = 1 - 2*s/(1 - ρ). Anybody knows how to compute the integral analytically?

The two parameters rho and s are real numbers between 0<s<1 && 0<rho <0.1 && 0<s+rho<1.

f = -((1 - 2*s + s^2 - ρ)^2/(-1 + s)^2)^(3/2) - ((1 - 2*s + s^2 - ρ)^2/(-1 + s)^2)^(3/2)*ρ + (12*z*ρ*(1 - 2*s + s^2 + ρ))/(-1 + s) - 
 (3*Sqrt[(1 - 2*s + s^2 - ρ)^2/(-1 + s)^2]*ρ*(1 - 2*s + s^2 + ρ))/(-1 + s) + 
 (3*z^2*Sqrt[(1 - 2*s + s^2 - ρ)^2/(-1 + s)^2]*ρ*(1 - 2*s + s^2 + ρ))/(-1 + s) + (12*z*ρ^2*(1 - 2*s + s^2 + ρ))/(-1 + s) + 
 (6*z*ρ*(1 - 2*s + s^2 + ρ)^2)/(-1 + s)^2 - (z*(1 - 2*s + s^2 + ρ)^3)/(-1 + s)^3 - (z*ρ*(1 - 2*s + s^2 + ρ)^3)/(-1 + s)^3 - 
 (3*z*(1 - 2*s + s^2 + ρ)^4)/(4*(-1 + s)^4) - (3*Sqrt[(1 - 2*s + s^2 - ρ)^2/(-1 + s)^2]*(1 - 2*s + s^2 + ρ)*
   (1 - 4*s + 6*s^2 - 4*s^3 + s^4 + ρ^2))/(8*(-1 + s)^3) - 
 (3*z^2*Sqrt[(1 - 2*s + s^2 - ρ)^2/(-1 + s)^2]*(1 - 2*s + s^2 + ρ)*(1 - 4*s + 6*s^2 - 4*s^3 + s^4 + ρ^2))/(8*(-1 + s)^3) + 
 (24*z*ρ*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ]))/(-1 + z^2) + (24*z*ρ^2*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ]))/
  (-1 + z^2) - (24*z*ρ*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2)/(-1 + z^2)^2 - 
 (8*z*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^3)/(-1 + z^2)^3 - (8*z*ρ*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^3)/
  (-1 + z^2)^3 + (12*z*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^4)/(-1 + z^2)^4 - 
 (12*ρ*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])*Sqrt[-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2])/
  (-1 + z^2) + (12*z^2*ρ*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])*
   Sqrt[-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2])/(-1 + z^2) + 
 8*(-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2)^(3/2) + 
 8*ρ*(-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2)^(3/2) - 
 (3*(1 - s - z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])*Sqrt[-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2]*
   (-ρ + (2*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2)/(-1 + z^2)^2))/(1 - z^2) - 
 (3*z^2*(1 - s - z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])*Sqrt[-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2]*
   (-ρ + (2*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2)/(-1 + z^2)^2))/(1 - z^2) - 
 15*ρ^2*ArcTanh[(-1 + 2*s - s^2 - ρ)/((-1 + s)*Sqrt[(1 - 2*s + s^2 - ρ)^2/(-1 + s)^2])] + 
 9*z^2*ρ^2*ArcTanh[(-1 + 2*s - s^2 - ρ)/((-1 + s)*Sqrt[(1 - 2*s + s^2 - ρ)^2/(-1 + s)^2])] + 
 15*ρ^2*ArcTanh[(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])/
    ((-1 + z^2)*Sqrt[-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2])] - 
 9*z^2*ρ^2*ArcTanh[(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])/
    ((-1 + z^2)*Sqrt[-ρ + (-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*ρ])^2/(-1 + z^2)^2])];

An example for the numerical integration would be:

zcut = 1 - 2*s/(1 - ρ);
NIntegrate[
 f*UnitStep[z - zcut]/. ρ -> 1.15^2/4.6^2 /. 
  s -> 0.4, {z, 0, 1}]

Result: 0.0714195

Answer 1

This is also an extended comment. I am hoping that the author of the OP will find it useful.

Note also, that the approach taken here is similar in spirit to the previous one by @lxndr . Since Integrate is listable and Mathematica refuses for her own reasons to do the full thing at once, let's try to help her as much as possible.

I am using f directly taken from the OP. Also the lower cut-off is given by

zcut = 1 - 2*s/(1 - ρ);

Ι am sitting on $Version

"13.0.0 for Mac OS X ARM (64-bit) (December 3, 2021)"

---

1. We split the integrand

   flist = List @@ Expand[f];

and we check how many terms we got

flist // Length

270

Nice!

Now, we invoke our friend Rubi

Get["Rubi`"]

and for completeness

$RubiVersion

"Rubi 4.16.1.0"

Before we proceed we break down the general idea. From those 270 terms in the expanded version of the integrand, Mathematica does some of them without any complaints, while she finds it difficult to deal with some others. For instance for the 21$^{st}$ term she complains

Integrate[flist[[21]], {z, zcut, 1}] // FullSimplify

Her issues are close to 1. So, we can subtract something from said value, say 1-a with $a>0$, integrate to the new point, and do the Series expansion as $a \rightarrow 0$. In other words,

Assuming[a > 0, Integrate[flist[[21]], {z, zcut, 1 - a}]] // 
   Series[#, {a, 0, 0}] & // Normal // FullSimplify

We can check the validity of this approach with a term for which she does not complain at all. Let's do the first term

Integrate[flist[[1]], {z, zcut, 1}] // FullSimplify
Assuming[a > 0, Integrate[flist[[1]], {z, zcut, 1 - a}]] // 
   Series[#, {a, 0, 0}] & // Normal // FullSimplify

both of the above return

Now, I am not sure if this approach agrees with the taste of the author of the OP, as it has not been addressed yet. I hope that they don't have any issues with this.

Edited bit 1:

How to go about it in general? Well, let's do an easy case. The first two terms. Integrate separately, do the Series around a goes to 0 and then sum the individual contributions.

For a list that looks like this

rubilist = 
 Table[Assuming[a > 0, Int[flist[[ii]], {z, zcut, 1 - a}]] // 
   FullSimplify, {ii, 1, 2}]

you can execute

FullSimplify[Total@Series[rubilist, {a, 0, 0}] // Normal]

that gives

End of edited bit 1.

---

WHY SHOULD WE USE RUBI?

Well, we do have 270 terms to consider. The following, hopefully, is a persuasive result.

listtestmma = 
 RepeatedTiming[
  Table[Assuming[a > 0, Integrate[flist[[ii]], {z, zcut, 1 - a}]] // 
     FullSimplify, {ii, 1, 21}];]

 listtestrubi = 
 RepeatedTiming[
  Table[Assuming[a > 0, Int[flist[[ii]], {z, zcut, 1 - a}]] // 
     FullSimplify, {ii, 1, 21}];]

So, it is clear how much better Rubi does in this particular example.

As a matter of fact, I checked the first 150 terms on my laptop

RepeatedTiming[
 Table[Assuming[a > 0, Int[flist[[ii]], {z, zcut, 1 - a}]] // 
    FullSimplify, {ii, 1, 150}];]
rubilist = 
  Table[Assuming[a > 0, Int[flist[[ii]], {z, zcut, 1 - a}]] // 
    FullSimplify, {ii, 1, 150}];

Practical note: if you follow this path, please take into consideration that Rubi does not fancy UnitStep as she does not have any rules for it, and that's why I re-wrote the integrand in the way I did.

I understand that this is not the end of the computation, but I am hoping that this is useful suggestion.

---

Edit 2: pushing Rubi and the suggested method or how I learned to stop worrying and took care of the first 177 terms.

This leaves 93 terms to be considered.

Note that the next term can also be evaluated. I just had to drop FullSimplify as it was taking forever. I used // TrigExpand // Simplify. Have to revisit.

rubilist1 = 
  Table[Assuming[a > 0, Int[flist[[ii]], {z, zcut, 1 - a}]] // 
    FullSimplify, {ii, 1, 177}];

Time to apply the proposed computational method to the above list

Simplify[
 PowerExpand[Factor[Total@Series[rubilist1, {a, 0, 0}] // Normal]]]

returns

For completeness, I am providing the results of the individual integrations

rubilist1 // TableForm

Edit 3: FullSimplify is Listable and maybe we can attempt to do so at the level of the integrand. 206 out of the 270 terms down. 64 to go

FullSimplify was mentioned but never addressed. It is not that difficult to be performed at the level of the integrand to try and simplify as much as possible before feeding something into Integrate

RepeatedTiming[
 TrigToExp@
   FullSimplify[Factor[flist[[1 ;; Length@flist]]], 
    0 < s < 1 && 0 < \[Rho] < 0.1];]

{1.04319, Null}

fsflist = 
  TrigToExp@
   FullSimplify[Factor[flist[[1 ;; Length@flist]]], 
    0 < s < 1 && 0 < \[Rho] < 0.1];

Then we can sum all the elements as

Total@fsflist

and hope that this integral will finish. However, this took more time than I could spare both using Rubi and without it - I have some upcoming deadlines, sorry.

So, this takes me back to the cheeky approach of computing the integration of each individual term in the list where we performed FullSimplify

In other words we perform

rubilistfs = 
  Table[Assuming[a > 0, Int[fsflist[[ii]], {z, zcut, 1 - a}]], {ii, 1,
     206}];

it takes a bit but it finishes. And this is as far as I have checked.

Answer 2

Just a hint, not a full solution: It seems to me that you need to simplify f so that some non converging terms cancel.

Try flist = List @@ Expand[f];.

Entry flist[[17]] is the first one that does not integrate, its value is (12 z)/(-1 + z^2)^4. For z->1 this gets infinite. I think that's the reason for Integrate[] to fail.

Since numerical integration converges there should be a counterterm canceling this. Simplify[] or FullSimplify[] should find and eliminate these, at least if the right assumptions are given. Unfortunately this simplification seems to be very time consuming - I stopped after some time of calculation.

Answer 3

To analyze function $F(\rho,s)\int _{zcut}^1fdz$ we need to know is this function analytical on $0<\rho <1$ and $0<s<1$. For this we map interval $(zcut,1)$ to unit interval with using variable $y=(z-zcut)/(1-zcut)$. Let us define function g=f/.{z -> y (1 - zcut) + zcut, s -> (1 - zcut) (1 - \[Rho])/2} then function F can be defined in numerical form as follows

FF[x1_?NumericQ, 
  x2_?NumericQ] := (1 - x2) NIntegrate[
   g /. {\[Rho] -> x1, zcut -> x2}, {y, 0, 1}]

Let check that FF is limited at the origin

FF[$MachineEpsilon, $MachineEpsilon]

Out[]= 0.132812 

With this definition we can plot function FF

Plot3D[FF[x1, x2], {x1, $MachineEpsilon, 1/2}, {x2, $MachineEpsilon, 
  1/2}, AxesLabel -> {"\[Rho]", "zcut", ""}, PlotTheme -> "Marketing",
  MeshStyle -> White, ColorFunction -> "Rainbow"]

Function F looks like analytical around origin, and therefore we have a chance to get some expression to describe F. To do this we expand function g in series and define several functions

ff = Normal[Series[g, {zcut, 0, 3}]];
coef=CoefficientList[ff, zcut];
F0[x1_?NumericQ] := NIntegrate[coef[[1]] /. \[Rho] -> x1, {y, 0, 1}]

data0 = Table[{x1, F0[x1]}, {x1, $MachineEpsilon, 1/2, .01}];

nlm0 = NonlinearModelFit[data0, 
  a + b x + c x^2 + d x^3 + e x^4 + h x^5, {a, b, c, d, e, h}, x];
 F1[x1_?NumericQ] := 
 NIntegrate[coef[[2]] /. \[Rho] -> x1, {y, 0, 1}, AccuracyGoal -> 5, 
   PrecisionGoal -> 4] // Quiet

data1 = Table[{x1, F1[x1]}, {x1, $MachineEpsilon, 1/2, .01}];
nlm1 = NonlinearModelFit[data1, 
  a + b x + c x^2 + d x^3 + e x^4 + h x^5, {a, b, c, d, e, h}, x];

F2[x1_?NumericQ] := NIntegrate[coef[[3]] /. \[Rho] -> x1, {y, 0, 1}];
    
data2 = Table[{x1, F2[x1]}, {x1, $MachineEpsilon, 1/2, .01}];

nlm2 = NonlinearModelFit[data2, 
  a + b x + c x^2 + d x^3 + e x^4 + h x^5, {a, b, c, d, e, h}, x];

F3[x1_?NumericQ] := 
 NIntegrate[coef[[4]] /. \[Rho] -> x1, {y, 0, 1}, AccuracyGoal -> 5, 
   PrecisionGoal -> 4] // Quiet

data3 = Table[{x1, F3[x1]}, {x1, $MachineEpsilon, 1/2, .01}];
nlm3 = NonlinearModelFit[data3, 
  a + b x + c x^2 + d x^3 + e x^4 + h x^5, {a, b, c, d, e, h}, x];

Finally we define analytical function described F with using series

fan[x_, x2_] := (1 - x2) (nlm0[x] + x2 nlm1[x] + x2^2 nlm2[x] + 
    x2^3 nlm3[x])

Visualization

Plot3D[fan[x1, x2], {x1, 0, 1/2}, {x2, 0, 
  1/2}, AxesLabel -> {"\[Rho]", "zcut", ""}, PlotTheme -> "Marketing",
  MeshStyle -> White, ColorFunction -> "Rainbow"]

Note, that fan looks very similar to function FF shown in Figure 1. We can express it in variables $\rho ,s$ and plot with using restrictions $\rho+s<1$

fff=((1 - x2) (nlm0[x] + x2 nlm1[x] + x2^2 nlm2[x] + 
      x2^3 nlm3[x]) // Normal) /. {x2 -> 1 - 2*s/(1 - \[Rho]), 
  x -> \[Rho]}



(* Out[]= 1/(
     1 - \[Rho]) 2 s (0.132915 - 0.85593 \[Rho] + 2.10123 \[Rho]^2 - 
       2.01754 \[Rho]^3 - 0.195384 \[Rho]^4 + 
       1.13491 \[Rho]^5 + (1 - (2 s)/(1 - \[Rho]))^2 (-0.1946 + 
          0.0518435 \[Rho] + 5.52769 \[Rho]^2 - 20.5944 \[Rho]^3 + 
          29.2689 \[Rho]^4 - 15.0774 \[Rho]^5) + (1 - (2 s)/(
          1 - \[Rho])) (-0.0472983 + 1.05398 \[Rho] - 4.94447 \[Rho]^2 + 
          9.34393 \[Rho]^3 - 7.08349 \[Rho]^4 + 1.14942 \[Rho]^5) + (1 - (
          2 s)/(1 - \[Rho]))^3 (0.0307526 + 0.112229 \[Rho] - 
          6.21303 \[Rho]^2 + 30.6017 \[Rho]^3 - 55.892 \[Rho]^4 + 
          35.9643 \[Rho]^5))*)


Plot3D[ If[\[Rho] + s < 1, fff, Nothing], {\[Rho], 0, 1/2}, {s, 0, 1},
  AxesLabel -> {"\[Rho]", "s", ""}, PlotTheme -> "Marketing", 
 MeshStyle -> White, ColorFunction -> "Rainbow", PlotRange -> All]

Update 1. We can improve code with using more series coefficient and make Do loop

g = -((1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2)^(3/
        2) - ((1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2)^(3/
        2)*\[Rho] + (12*z*\[Rho]*(1 - 2*s + s^2 + \[Rho]))/(-1 + 
       s) - (3*Sqrt[(1 - 2*s + s^2 - \[Rho])^2/(-1 + 
            s)^2]*\[Rho]*(1 - 2*s + s^2 + \[Rho]))/(-1 + s) + (3*z^2*
       Sqrt[(1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2]*\[Rho]*(1 - 2*s + 
         s^2 + \[Rho]))/(-1 + s) + (12*
       z*\[Rho]^2*(1 - 2*s + s^2 + \[Rho]))/(-1 + s) + (6*
       z*\[Rho]*(1 - 2*s + s^2 + \[Rho])^2)/(-1 + 
        s)^2 - (z*(1 - 2*s + s^2 + \[Rho])^3)/(-1 + 
        s)^3 - (z*\[Rho]*(1 - 2*s + s^2 + \[Rho])^3)/(-1 + s)^3 - (3*
       z*(1 - 2*s + s^2 + \[Rho])^4)/(4*(-1 + s)^4) - (3*
       Sqrt[(1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2]*(1 - 2*s + 
         s^2 + \[Rho])*(1 - 4*s + 6*s^2 - 4*s^3 + 
         s^4 + \[Rho]^2))/(8*(-1 + s)^3) - (3*z^2*
       Sqrt[(1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2]*(1 - 2*s + 
         s^2 + \[Rho])*(1 - 4*s + 6*s^2 - 4*s^3 + 
         s^4 + \[Rho]^2))/(8*(-1 + s)^3) + (24*
       z*\[Rho]*(-1 + s + 
         z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]]))/(-1 + z^2) + (24*
       z*\[Rho]^2*(-1 + s + 
         z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]]))/(-1 + z^2) - (24*
       z*\[Rho]*(-1 + s + 
          z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2)/(-1 + 
        z^2)^2 - (8*
       z*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^3)/(-1 +
         z^2)^3 - (8*
       z*\[Rho]*(-1 + s + 
          z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^3)/(-1 + 
        z^2)^3 + (12*
       z*(-1 + s + z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^4)/(-1 +
         z^2)^4 - (12*\[Rho]*(-1 + s + 
         z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])*
       Sqrt[-\[Rho] + (-1 + s + 
             z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
             z^2)^2])/(-1 + z^2) + (12*
       z^2*\[Rho]*(-1 + s + 
         z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])*
       Sqrt[-\[Rho] + (-1 + s + 
             z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
             z^2)^2])/(-1 + z^2) + 
    8*(-\[Rho] + (-1 + s + 
            z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
            z^2)^2)^(3/2) + 
    8*\[Rho]*(-\[Rho] + (-1 + s + 
            z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
            z^2)^2)^(3/
        2) - (3*(1 - s - z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])*
       Sqrt[-\[Rho] + (-1 + s + 
             z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
             z^2)^2]*(-\[Rho] + (2*(-1 + s + 
               z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2)/(-1 + 
             z^2)^2))/(1 - z^2) - (3*
       z^2*(1 - s - z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])*
       Sqrt[-\[Rho] + (-1 + s + 
             z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
             z^2)^2]*(-\[Rho] + (2*(-1 + s + 
               z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2)/(-1 + 
             z^2)^2))/(1 - z^2) - 
    15*\[Rho]^2*
     ArcTanh[(-1 + 2*s - s^2 - \[Rho])/((-1 + s)*
         Sqrt[(1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2])] + 
    9*z^2*\[Rho]^2*
     ArcTanh[(-1 + 2*s - s^2 - \[Rho])/((-1 + s)*
         Sqrt[(1 - 2*s + s^2 - \[Rho])^2/(-1 + s)^2])] + 
    15*\[Rho]^2*
     ArcTanh[(-1 + s + 
         z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])/((-1 + z^2)*
         Sqrt[-\[Rho] + (-1 + s + 
               z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
               z^2)^2])] - 
    9*z^2*\[Rho]^2*
     ArcTanh[(-1 + s + 
         z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])/((-1 + z^2)*
         Sqrt[-\[Rho] + (-1 + s + 
               z*Sqrt[1 - 2*s + s^2 + (-1 + z^2)*\[Rho]])^2/(-1 + 
               z^2)^2])] /. {z -> y (1 - zcut) + zcut, 
    s -> (1 - zcut) (1 - \[Rho])/2};

ff = Normal[Series[g, {zcut, 0, 4}]];

coef = CoefficientList[ff, zcut]; nn = Length[coef];

Do[data[i] = 
  Table[{x1, 
    NIntegrate[coef[[i]] /. \[Rho] -> x1, {y, 0, 1}, 
      AccuracyGoal -> 5, PrecisionGoal -> 4] // 
     Quiet}, {x1, $MachineEpsilon, 1/2, .01}]; 
 nlm[i] = NonlinearModelFit[data[i], 
   a + b x + c x^2 + d x^3 + e x^4 + h x^5, {a, b, c, d, e, h}, 
   x];, {i, nn}]

fan[x_, x2_] := (1 - x2) Sum[nlm[n][x] x2^(n - 1), {n, 1, nn}]
fff = ((1 - x2) Sum[nlm[n][x] x2^(n - 1), {n, 1, nn}] // 
    Normal);

fff //. {x2 -> 1 - 2*s/(1 - \[Rho]), x -> \[Rho]}

(*Out[]= 1/(
 1 - \[Rho]) 2 s (0.132915 - 0.85593 \[Rho] + 2.10123 \[Rho]^2 - 
   2.01754 \[Rho]^3 - 0.195384 \[Rho]^4 + 
   1.13491 \[Rho]^5 + (1 - (2 s)/(1 - \[Rho]))^4 (0.0542876 - 
      0.331013 \[Rho] + 6.4239 \[Rho]^2 - 35.6162 \[Rho]^3 + 
      73.6482 \[Rho]^4 - 52.5108 \[Rho]^5) + (1 - (2 s)/(
      1 - \[Rho]))^2 (-0.1946 + 0.0518435 \[Rho] + 5.52769 \[Rho]^2 - 
      20.5944 \[Rho]^3 + 29.2689 \[Rho]^4 - 15.0774 \[Rho]^5) + (1 - (
      2 s)/(1 - \[Rho])) (-0.0472983 + 1.05398 \[Rho] - 
      4.94447 \[Rho]^2 + 9.34393 \[Rho]^3 - 7.08349 \[Rho]^4 + 
      1.14942 \[Rho]^5) + (1 - (2 s)/(1 - \[Rho]))^3 (0.0307526 + 
      0.112229 \[Rho] - 6.21303 \[Rho]^2 + 30.6017 \[Rho]^3 - 
      55.892 \[Rho]^4 + 35.9643 \[Rho]^5))*)

We can plot final result to observe improvement around line $\rho +s=1$

Plot3D[ If[\[Rho] + s < 1, 
  fff //. {x2 -> 1 - 2*s/(1 - \[Rho]), x -> \[Rho]}, 
  Nothing], {\[Rho], 0, 1/2}, {s, 0, 1}, 
 AxesLabel -> {"\[Rho]", "s", ""}, PlotTheme -> "Marketing", 
 MeshStyle -> White, ColorFunction -> "Rainbow", PlotRange -> All]

But direct computation with FF shows different picture

Plot3D[ If[x + s < 1, FF[x, 1 - 2*s/(1 - x)], 
  Nothing], {x, $MachineEpsilon, 1/2}, {s, $MachineEpsilon, 1}, 
 AxesLabel -> {"\[Rho]", "s", ""}, PlotTheme -> "Marketing", 
 MeshStyle -> White, ColorFunction -> "Rainbow", PlotRange -> All]

It means that there is different criterium 1 - 2*s + s^2 - \[Rho] >= 0 follows from g definition. With this criterium we have picture very similar as above

Plot3D[ If[1 - 2*s + s^2 - \[Rho] >= 0, 
  fff //. {x2 -> 1 - 2*s/(1 - \[Rho]), x -> \[Rho]}, 
  Nothing], {\[Rho], 0, 1/2}, {s, 0, 1}, 
 AxesLabel -> {"\[Rho]", "s", ""}, PlotTheme -> "Marketing", 
 MeshStyle -> White, ColorFunction -> "Rainbow", PlotRange -> All, 
 PlotPoints -> 150]