Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- from scipy import integrate
- from scipy.special import jacobi
- from tabulate import tabulate
- a = -1.0
- b = 1.0
- Q1 = np.poly1d([1, 1.5, 0, -14.5])
- Q2 = np.poly1d([1, 2, -11])
- def f(x):
- return 1. - x / 2.
- def points(n):
- return np.linspace(a, b, n)
- def colloc_points(n):
- return [np.cos(((2 * i - 1) * np.pi) / (2 * n)) for i in range(1, n + 1)]
- def w(n):
- if n == 1:
- return [Q1]
- elif n == 2:
- return [Q1, Q2]
- else:
- jacobi_tail = [np.poly1d([1, 0, -2, 0, 1]) *
- jacobi(i - 3, 2, 2) for i in range(3, n + 1)]
- return [Q1, Q2] + jacobi_tail
- def L_diff(w, x):
- w_2 = np.polyder(w, 2)
- w_1 = np.polyder(w, 1)
- ans = -1 / (3 + np.sin(x)) * w_2(x) - x * w_1(x) + (1 + x) * w(x)
- return ans
- def L_colloc(n):
- L = np.zeros((n, n))
- g = np.zeros(n)
- W = w(n)
- X = colloc_points(n)
- for i in range(0, n):
- for j in range(0, n):
- L[i, j] = L_diff(W[j], X[i])
- g[i] = f(X[i])
- return L, g
- def L_gal(n):
- L = np.zeros((n, n))
- g = np.zeros(n)
- W = w(n)
- for j in range(0, n):
- for i in range(0, n):
- L[j, i] = integrate.quad(lambda x: W[j](x) * L_diff(W[i], x), a, b)[0]
- g[j] = integrate.quad(lambda x: W[j](x) * f(x), a, b)[0]
- return L, g
- def coeff(l, n):
- L, g = l(n)
- return np.linalg.solve(L, g)
- def solution(W, C, n):
- return lambda x: np.dot(np.array([W[i](x) for i in range(0, n)]), C)
- def print_table(n, a, b, l):
- print((' ' * 6) + "|" + (' ' * 14) + 'y_n(x)' + (' ' * 16))
- table = []
- for i in range(3, n + 1):
- W = w(i)
- C = coeff(l, i)
- sol = solution(W, C, i)
- table.append([i, sol(a), sol((a + b) / 2), sol(b)])
- head = ['x = {0}'.format(a), 'x = {0}'.format((a + b) / 2), 'x = {0}'.format(b)]
- print(tabulate(table, headers=['n'] + head * 2, tablefmt="grid", floatfmt="f"))
- print('Galerkin method')
- print('')
- print_table(7, a, b, L_gal)
- print('')
- print('Collocation method')
- print('')
- print_table(7, a, b, L_colloc)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement