#!/usr/bin/env python3 """Chapter 10, code 2: a linear BVP with different boundary conditions. Solve y'' + p(x)*y' + q(x)*y = r(x) on a uniform grid. Method: centered finite differences and the Thomas algorithm. Derivative boundaries use ghost-point elimination; endpoints are unknowns. Example: dimensionless heat conduction in a rod, T'' = -sin(pi*x). Default: T(0) = 0, T'(1) = 0 (an insulated right end). y denotes T. Exact solutions are supplied for all four boundary keys. Requires NumPy and Matplotlib. Run: python3 ch10_boundary_conditions.py SciPy alternatives: solve_banded for a banded system, or spsolve for a sparse matrix. Include the selected boundary rows in the assembled system. To solve another problem, edit the inputs and coefficients below; replace exact_solution or set check_exact = False. """ import numpy as np import matplotlib.pyplot as plt def main(): # ------------------- Problem and grid: edit here ------------------ a = 0.0 b = 1.0 n = 20 key = 2 left_value = 0.0 right_value = 0.0 robin_alpha = 1.0 robin_beta = 1.0 check_exact = True make_plot = True # Boundary key and meaning of the two boundary data: # 1: y(a) = left_value; y(b) = right_value # 2: y(a) = left_value; y'(b) = right_value # 3: y'(a) = left_value; y(b) = right_value # 4: y(a) = left_value; # robin_alpha*y(b) + robin_beta*y'(b) = right_value # Derivatives mean dy/dx, not the outward normal derivative. # For a rod's convecting right end, key 4 expresses # k*T'(b) + H*T(b) = H*T_environment. # Use dimensionless conductivity k and heat-transfer coefficient H. # ----------------------------------------------------------------- if n < 2 or n != int(n) or b <= a: raise ValueError('Use integer n >= 2 and b > a.') if key not in (1, 2, 3, 4): raise ValueError('Choose a boundary key from 1 to 4.') if key == 4: if robin_beta == 0: raise ValueError('Use key 1 for a fixed right-end value.') scale = max(abs(robin_alpha*(b-a)), abs(robin_beta)) if check_exact and abs(robin_alpha*(b-a)+robin_beta) <= np.finfo(float).eps*scale: raise ValueError('The supplied example has no unique solution for these Robin data.') h = (b-a)/n x = np.linspace(a, b, n+1) y = finite_difference(coefficients, n, x, key, left_value, right_value, robin_alpha, robin_beta) print('Linear boundary value problem: boundary-condition choices') print(f'Boundary key: {key:8d}') print(f'Grid intervals: {n:8d}') print(f'Grid spacing: {h:12.6f}') descriptions = { 1: 'Fixed values at both ends.', 2: 'Fixed left value; prescribed right derivative.', 3: 'Prescribed left derivative; fixed right value.', 4: 'Fixed left value; right Robin condition.'} print(descriptions[key]+'\n') filename = 'ch10_boundary_conditions.dat' if check_exact: exact = exact_solution(x, a, b, key, left_value, right_value, robin_alpha, robin_beta) abs_error = np.abs(y-exact) with open(filename, 'w', encoding='utf-8') as data_file: if check_exact: print(f"{'x':>10} {'numerical y':>16} {'exact y':>16} {'abs error':>12}") data_file.write('# x numerical_y exact_y abs_error\n') else: print(f"{'x':>10} {'numerical y':>16}") data_file.write('# x numerical_y\n') for i in range(n+1): if check_exact: print(f'{x[i]:10.6f} {y[i]:16.8e} {exact[i]:16.8e} {abs_error[i]:12.4e}') row = (x[i], y[i], exact[i], abs_error[i]) else: print(f'{x[i]:10.6f} {y[i]:16.8e}') row = (x[i], y[i]) data_file.write(' '.join(f'{value:.16e}' for value in row)+'\n') print() if check_exact: print(f'Largest grid error: {np.max(abs_error):12.4e}') print(f'Data file: {filename}') if make_plot: panels = 2 if check_exact else 1 fig, axes = plt.subplots(panels, 1, figsize=(7, 7 if check_exact else 4), squeeze=False) ax = axes[0, 0] ax.plot(x, y, 'ko-', linewidth=1, markersize=4, label='Finite differences') if check_exact: x_plot = np.linspace(a, b, 401) ax.plot(x_plot, exact_solution(x_plot, a, b, key, left_value, right_value, robin_alpha, robin_beta), 'k--', linewidth=1.2, label='Exact solution') ax.legend() ax.set_xlabel('x') ax.set_ylabel('y(x)') ax.set_title('Linear boundary value problem') ax.grid(True) if check_exact: ax = axes[1, 0] ax.plot(x, abs_error, 'ko-', linewidth=1, markersize=4) ax.set_xlabel('x') ax.set_ylabel('Absolute error') ax.grid(True) fig.tight_layout() plt.show() def coefficients(x): """Supply p, q, r in y'' + p(x)*y' + q(x)*y = r(x). Here y is dimensionless temperature in a uniformly conducting rod. """ p = 0.0 q = 0.0 r = -np.sin(np.pi*x) return p, q, r def exact_solution(x, a, b, key, left_value, right_value, alpha, beta): """For y'' = -sin(pi*x), y = sin(pi*x)/pi^2 + c*(x-a) + d. Determine c and d from the selected boundary conditions. Default key 2 on [0,1]: y = sin(pi*x)/pi^2 + x/pi. """ pa = np.sin(np.pi*a)/np.pi**2 pb = np.sin(np.pi*b)/np.pi**2 if key == 1: d = left_value-pa c = (right_value-pb-d)/(b-a) elif key == 2: d = left_value-pa c = right_value-np.cos(np.pi*b)/np.pi elif key == 3: c = left_value-np.cos(np.pi*a)/np.pi d = right_value-pb-c*(b-a) else: # key 4 d = left_value-pa c = (right_value-alpha*(pb+d)-beta*np.cos(np.pi*b)/np.pi)/(alpha*(b-a)+beta) return np.sin(np.pi*x)/np.pi**2+c*(x-a)+d def finite_difference(fcn, n, x, key, left_value, right_value, alpha, beta): """Assemble all n+1 unknowns, including the endpoints. Interior equations are identical to Code 1. A fixed endpoint gets a row containing y = prescribed value. For a derivative boundary, introduce a ghost point outside the grid: y_ghost_left = y_next - 2*h*y'(a) y_ghost_right = y_previous + 2*h*y'(b). Substitute into the centered ODE at the endpoint to remove the ghost. For the right Robin condition, use y'(b) = (right_value-alpha*y(b))/beta. Python index 0 is the left endpoint; index n is the right endpoint. """ h = (x[-1]-x[0])/n lower = np.zeros(n+1) diagonal = np.zeros(n+1) upper = np.zeros(n+1) rhs = np.zeros(n+1) for i in range(1, n): p, q, r = fcn(x[i]) lower[i] = 1-h*p/2 diagonal[i] = -2+h*h*q upper[i] = 1+h*p/2 rhs[i] = h*h*r # Left boundary: derivative for key 3; fixed value otherwise. if key == 3: p, q, r = fcn(x[0]) diagonal[0] = -2+h*h*q upper[0] = 2 rhs[0] = h*h*r+(2*h-h*h*p)*left_value else: diagonal[0] = 1 rhs[0] = left_value # Right boundary: fixed value, derivative, or Robin condition. if key in (1, 3): diagonal[n] = 1 rhs[n] = right_value elif key == 2: p, q, r = fcn(x[n]) lower[n] = 2 diagonal[n] = -2+h*h*q rhs[n] = h*h*r-(2*h+h*h*p)*right_value else: # key 4 p, q, r = fcn(x[n]) factor = (2*h+h*h*p)/beta lower[n] = 2 diagonal[n] = -2+h*h*q-factor*alpha rhs[n] = h*h*r-factor*right_value return thomas(lower, diagonal, upper, rhs) def thomas(lower, diagonal, upper, rhs): """Forward elimination and backward substitution for a tridiagonal system. Copies preserve the caller's diagonal and rhs arrays. No pivoting is used. A zero or near-zero pivot stops the calculation; this may require a solver with pivoting or indicate a singular BVP. """ diagonal = diagonal.copy() rhs = rhs.copy() m = len(rhs) y = np.zeros(m) for i in range(m): scale = max(abs(lower[i]), abs(diagonal[i]), abs(upper[i])) if abs(diagonal[i]) <= np.finfo(float).eps*scale: raise RuntimeError('Thomas algorithm: zero or near-zero pivot.') if i == m-1: break multiplier = lower[i+1]/diagonal[i] diagonal[i+1] -= multiplier*upper[i] rhs[i+1] -= multiplier*rhs[i] y[-1] = rhs[-1]/diagonal[-1] for i in range(m-2, -1, -1): y[i] = (rhs[i]-upper[i]*y[i+1])/diagonal[i] return y if __name__ == '__main__': main()