Gilbert Strang explains the matrix square root of the second difference matrix [here](http://math.mit.edu/highdegree/TplusH.pdf).

In particular,

<pre>
import numpy
import scipy.linalg

def f(N, p):
    a = 2 * (-1)**(p-1)
    b = 1/numpy.tan((numpy.pi * (2*p - 1)) / (4 + 4*N))
    c = 1/numpy.tan((numpy.pi * (2*p + 1)) / (4 + 4*N))
    return a - b + c

N = 10000
K_tri = numpy.eye(N) - numpy.eye(N, k=1)
K = K_tri + K_tri.T

s = numpy.arange(N)
T = scipy.linalg.toeplitz(f(N, s))
H = scipy.linalg.hankel(f(N, s+2), f(N, N+1+s))
K_sqrt = (0.5 / (N+1)) * (T - H)

print(numpy.max(numpy.abs(numpy.dot(K_sqrt, K_sqrt) - K)))
</pre>

<pre>
2.14717132963e-13
</pre>