Inconsistencies in time calculations for a Brachistochrone curve with friction

Question

I am trying to find the time taken by a bead to slide down the Brachistochrone curve (between 0,0 and 1,-1) under kinetic friction, but I'm getting inconsistent answers with analytical and numerical methods.

This article on Wolfram MathWorld solves the problem analytically, giving the equations:

$$ x = a[(\theta-\sin{\theta})+\mu(1-\cos{\theta})] $$

$$ y = a[(1-\cos{\theta})+\mu(\theta+\sin{\theta})] $$

Using the same method in this ProofWiki article with modified equations gives me this integral:

$$ \int_{0}^{2.84168} \sqrt{\frac{0.2837748856^2 \left( (1 - \cos(t) + 0.5 \sin(t))^2 + (\sin(t) + 0.5 (1 + \cos(t)))^2 \right)}{2 \times 9.8 \times 0.2837748856 \times (1 - \cos(t) + 0.5 (t + \sin(t)))}} dt = 0.586767 $$

(I'm apprehensive about this calculation; it seems too close (if that's a thing) to the frictionless time value of $0.583196$.)

However, when I run some optimization code to compute the Brachistochrone with friction, I consistently get values in the $0.66-0.67$ range. I got the code from a YouTube video and modified it to factor in friction.

Code Block 1:

g = 9.8
mu_k = 0.5

# Define functions for finding time between two points and total time for the path
def dtime_fric(x1, x2, y1, y2, v1):
    p1 = np.array([x1, y1])
    p2 = np.array([x2, y2])
    yp = y2 - y1
    xp = x2 - x1
    s = np.linalg.norm(p2 - p1)
    a = -g * yp / s
    acc_fric = mu_k * (g * xp / s)
    a -= acc_fric

    if yp == 0:
        tp = s / v1
        v2 = v1
    else:
        tp = (-v1 + np.sqrt(v1**2 + 2 * a * s)) / a
        v2 = v1 + a * tp
    return tp, v2

def calc_time_fric(xx, pp):
    tt = 0
    vc = 0
    for i in range(1, len(pp)):
        ttemp, vnext = dtime_fric(xx[i-1], xx[i], pp[i-1], pp[i], vc)
        tt += ttemp
        vc = vnext
    return tt

# Define initial points and parameters
p1x, p1y = 0, 0
p2x, p2y = 1, -1
N = 50
dx = (p2x - p1x) / N
dy = (p2y - p1y) / N

# Initialize arrays
p = []
pn = []
x = []
n = 0

# Make points
while n <= N:
    x.append(p1x + n * dx)
    p.append(p1y + n * dy)
    pn.append(p1y + n * dy)
    n += 1

# Initialize variables
nn = 0
NN = 100_000
t_old = calc_time_fric(x, p)

# Random optimization loop
while nn < NN:
  if nn%20000 == 0:
    print("Running iteration: ", nn)
  # for i in range(1, len(pn)-1):
  for i in range(len(pn)-2, 0, -1):
    pn[i] = p[i] + 0.002 * (2 * np.random.random() - 1)
  t_new = calc_time_fric(x, pn)
  if t_new < t_old:
    t_old = calc_time_fric(x, pn)
    for i in range(1, len(p)-1):
      p[i] = pn[i]
  nn += 1

I'm probably doing this wrong but when I use the Brach w/ friction equations and run a time calculation function from this Github repo on it, I get $0.5864$, which is eerily close to the value I got from that first integral.

Code Block 2:

def make_curve_cycloid_friction():
    N = 100_000
    a = a_fric
    t = t_fric
    list_t = np.linspace(0.0, t, N)
    list_x = a * (list_t - np.sin(list_t) + mu_k * (1 - np.cos(list_t)))
    list_y = a * (-1.0 + np.cos(list_t) - mu_k * (list_t + np.sin(list_t)))
    return list_x, list_y

def calc_time_friction(list_x, list_y, g):
    time = 0.0
    len_x = len(list_x)
    list_dx = np.array([list_x[i+1] - list_x[i] for i in range(len_x-1)])
    list_dy = np.array([list_y[i+1] - list_y[i] for i in range(len_x-1)])
    list_time = np.array([np.sqrt(((list_dx[i])**2 + (list_dy[i])**2) / (0.5*np.abs(list_y[i] + list_y[i+1]))) for i in range(len_x-1)])
    time = np.sum(list_time) / np.sqrt(2.0*g)
    return time

However, when I use the calc_time_fric from the code block 1 on the coordinate lists from make_curve_cycloid_friction, I get $0.6545$, which is close to the optimization curve time from code block 1.

So I'm getting two different "regions" of values from my calculations. Have I coded the friction logic wrong? Is my integral incorrect? Is one of my calc_time functions wrong? Or is the problem something else altogether?

EDIT 1: After having another look at the MathWorld Brach w/ friction derivation, I noticed that the formula for speed when there's friction is actually $$v=\sqrt{2g(y-\mu x)}$$ I changed the calc_time_friction function in code block 2 accordingly:

    list_time = np.array([np.sqrt(((list_dx[i])**2 + (list_dy[i])**2) / ( (0.5*np.abs(list_y[i] + list_y[i+1])) - mu_k * (0.5*np.abs(list_x[i] + list_x[i+1])))) for i in range(len_x-1)])
    time = np.sum(list_time) / np.sqrt(2.0*g)

and the value I am now getting from this calculation is $0.6541$. I think that this, along with the values from my optimization code and the calc_time_fric function from code block 1, suggests that the values in the $0.65-0.67$ range are correct, and there's something incorrect in the first integral.

EDIT 2: Alright I fixed the equation of $v$ in the integral and got $0.65447$.

$$ \int_{0}^{2.84168} \sqrt{\frac{0.2837748856^2 \left( (1 - \cos(t) + 0.5 \sin(t))^2 + (\sin(t) + 0.5 (1 + \cos(t)))^2 \right)}{2 \times 9.8 \times \left(0.2837748856 \left(1 - \cos(t) + 0.5 (t + \sin(t))\right) - 0.5 \left(0.2837748856 \left(t - \sin(t) + 0.5 (1 - \cos(t))\right)\right)\right)}} dt = 0.65447 $$

This basically resolves my question, I guess.