import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma, lpmv
nxs = 301
x = np.linspace(-1, 1, nxs)
nterms = [10, 50, 100]
u = np.zeros((len(nterms), nxs))
ntermsmax = nterms[-1]
phi = np.zeros((ntermsmax+1, nxs))
for m in range(ntermsmax+1):
P = lpmv(0, m, x)
phi[m, :] = P * np.sqrt((2 * m + 1) / 2)
flist = (np.sign(x) + 1) / 2
for i in range(len(nterms)):
k = 0
un0 = np.sqrt(np.pi) / (2 * gamma(1 - k / 2) * gamma((3 + k) / 2)) * np.sqrt((2 * k + 1) / 2)
u[i, :] += un0 * phi[k, :]
for k in range(1, nterms[i] + 1, 2):
un0 = np.sqrt(np.pi) / (2 * gamma(1 - k / 2) * gamma((3 + k) / 2)) * np.sqrt((2 * k + 1) / 2)
u[i, :] += un0 * phi[k, :]
plt.plot(x, flist, '-', x, u[0, :], '--', x, u[1, :], ':', x, u[2, :], '-.')
plt.legend(['H(x)', 'n=10', 'n=50', 'n=100'])
plt.show(block=False)
table = np.column_stack( (x, flist, u.T) )
with open("LegendreStep.dat", "w") as f:
np.savetxt(f, table, fmt='%f', header="table")