Figure 2.31:

Dependence of \left | c(j) \right | on j for the Legendre-Galerkin approximation of \eqref {eq:MWRexample} with n=10.

Code for Figure 2.31

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import numpy as np
import matplotlib.pyplot as plt

N = 10
A = np.zeros((N, N))
b = np.zeros((N, 1))

for i in range(N-2):
    for j in range(N):
        A[i, j] = 2. / (2 * j + 1) * (i == j)
        for k in range(j+1):
            if (j-k) % 2 == 0:
                A[i, j] += 4 * (j - k) * (j + k + 1) * (i == k)

for j in range(N):
    A[N-2, j] = (-1)**j
    A[N-1, j] = 1

for i in range(N-2):
    b[i] = -(i == 0) - 1 / 3 * (i == 1)

b[N-2] = 0
b[N-1] = 0
c = np.linalg.solve(A, b)

nxs = 1000
x = np.linspace(-1, 1, nxs)
soln = np.zeros(nxs)
phi = np.zeros((N, nxs))

for m in range(N):
    P = np.polynomial.Legendre.basis(m)
    phi[m, :] = P(x)

for j in range(N):
    for i in range(nxs):
        soln[i] += c[j] * phi[j, i]

Nex = 100
hex = 1 / (Nex + 1)
xex = np.linspace(0, 1, Nex+2)
uex = -xex + np.sin(xex) / np.sin(1)

plt.figure()
plt.plot((x + 1) / 2, soln, label='Numerical Solution')
plt.plot(xex, uex, 'r-', label='Exact Solution')
plt.legend()
plt.title('Numerical and Exact Solution')

ind = np.arange(1, N+1)
plt.figure()
plt.semilogy(ind, abs(c), 'o')
plt.title('Semilog Plot of Coefficients')

plt.show(block=False)

t = np.column_stack((ind, abs(c)))

with open("LegendreGalerkin1.dat", "w") as f:
    np.savetxt(f, t, fmt='%e', header="Legendre coefficients")