Return to Answer

Gilbert Strang explains the matrix square root of the second difference matrix here.

In particular,

import numpy
import scipy.linalg

def f(p):
    return (4 / numpy.pi) / ((1 - 2*p) * (1 + 2*p))

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(s))
H = -scipy.linalg.hankel(f(s+2), f(N+1-s))
K_sqrt = T + H

print(numpy.max(numpy.abs(numpy.dot(K_sqrt, K_sqrt) - K)))

7.90678633678e-12

Gilbert Strang explains the matrix square root of the second difference matrix here.

In particular,

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)))

2.14717132963e-13

Gilbert Strang explains the matrix square root of the second difference matrix here.

Gilbert Strang explains the matrix square root of the second difference matrix here.

In particular,

import numpy
import scipy.linalg

def f(p):
    return (4 / numpy.pi) / ((1 - 2*p) * (1 + 2*p))

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(s))
H = -scipy.linalg.hankel(f(s+2), f(N+1-s))
K_sqrt = T + H

print(numpy.max(numpy.abs(numpy.dot(K_sqrt, K_sqrt) - K)))

7.90678633678e-12

Gilbert Strang explains the matrix square root of the second difference matrix here.

Gilbert Strang explains the matrix square root of the second difference matrix here.