#!/usr/bin/env python3 """Chapter 10, code 3: a nonlinear boundary value problem. Solve y'' = f(x,y,y'), with y(a) = alpha and y(b) = beta. Method: centered finite differences and Newton iteration. Each Newton correction is found with the Thomas algorithm. Example: y'' = y^2 - sin(pi*x)^2 - pi^2*sin(pi*x), on [0,1]. y(0) = y(1) = 0; exact solution y = sin(pi*x). The forcing was chosen so that the exact solution is known. Requires NumPy and Matplotlib. Run: python3 ch10_nonlinear_bvp.py SciPy alternatives: solve_banded can solve each Newton system, or solve_bvp can solve the BVP directly by collocation using a first-order system and an initial solution guess. To solve another problem: 1. Edit the interval, endpoint values, grid, and initial guess in main. 2. Edit equation below, including its two partial derivatives. 3. Disable check_exact or replace exact_solution. """ import numpy as np import matplotlib.pyplot as plt def main(): # ---------------- Problem and controls: edit here ----------------- a = 0.0 b = 1.0 alpha = 0.0 # y(a) beta = 0.0 # y(b) n = 20 # Number of grid intervals tolerance = 1.0e-10 # Maximum absolute discrete-equation residual max_iterations = 30 check_exact = True make_plot = True # ----------------------------------------------------------------- if n < 2 or n != int(n) or b <= a: raise ValueError('Use integer n >= 2 and b > a.') if tolerance <= 0 or max_iterations < 1 or max_iterations != int(max_iterations): raise ValueError('Use positive tolerance and an integer iteration limit.') h = (b-a)/n x = np.linspace(a, b, n+1) # Initial guess: a straight line between the endpoint values. # For another initial guess, change this expression. y = alpha+(beta-alpha)*(x-a)/(b-a) y[0] = alpha y[-1] = beta print('Nonlinear boundary value problem: finite differences and Newton iteration') print(f'Grid intervals: {n:8d}') print(f'Grid spacing: {h:12.6f}') print(f'Newton tolerance: {tolerance:12.4e}') y, iterations, residual_norm = newton_bvp( equation, n, x, y, tolerance, max_iterations) print() filename = 'ch10_nonlinear_bvp.dat' if check_exact: exact = exact_solution(x) 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() print(f'Newton iterations: {iterations:8d}') print(f'Final residual: {residual_norm:12.4e}') 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), 'k--', linewidth=1.2, label='Exact solution') ax.legend() ax.set_xlabel('x') ax.set_ylabel('y(x)') ax.set_title('Nonlinear 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 equation(x, y, v): """Supply f, df/dy and df/dv for y'' = f(x,y,v), where v = y'. Partial derivatives treat x, y and v as independent variables. If f does not depend on v, set dfdv = 0. """ s = np.sin(np.pi*x) f = y*y-s*s-np.pi*np.pi*s dfdy = 2*y dfdv = 0.0 return f, dfdy, dfdv def exact_solution(x): """Exact solution for the supplied equation and default boundary values.""" return np.sin(np.pi*x) def newton_bvp(fcn, n, x, y, tolerance, max_iterations): """Solve the discrete equations using full Newton steps. At each interior point: R = (y_left-2*y_center+y_right)/h^2 - f(x,y_center,v), v = (y_right-y_left)/(2*h). Solve J*correction = -R. Multiplying the Newton equations by h^2 gives lower = 1+h*(df/dv)/2, diagonal = -2-h^2*(df/dy), upper = 1-h*(df/dv)/2, rhs = -h^2*R. Interior unknown i lies at x[i+1]. Endpoint corrections are zero. Stop when max|R| <= tolerance. This measures convergence of the discrete equations; grid refinement checks discretization error. An unsuitable initial guess may cause failure or convergence to a different solution when several solutions exist. """ y = y.copy() h = (x[-1]-x[0])/n m = n-1 lower = np.zeros(m) diagonal = np.zeros(m) upper = np.zeros(m) rhs = np.zeros(m) change = 0.0 print(f"{'iteration':>10} {'max residual':>18} {'max correction':>18}") for iterations in range(max_iterations+1): residual_norm = 0.0 for i in range(m): j = i+1 v = (y[j+1]-y[j-1])/(2*h) f, dfdy, dfdv = fcn(x[j], y[j], v) if not np.all(np.isfinite([f, dfdy, dfdv])): raise RuntimeError('Nonfinite equation or partial derivative; check initial guess.') residual = (y[j-1]-2*y[j]+y[j+1])/(h*h)-f residual_norm = max(residual_norm, abs(residual)) lower[i] = 1+h*dfdv/2 diagonal[i] = -2-h*h*dfdy upper[i] = 1-h*dfdv/2 rhs[i] = -h*h*residual # Correction shown was applied in the preceding iteration. print(f'{iterations:10d} {residual_norm:18.10e} {change:18.10e}') if residual_norm <= tolerance: return y, iterations, residual_norm if iterations == max_iterations: break lower[0] = 0 upper[-1] = 0 correction = thomas(lower, diagonal, upper, rhs) change = np.max(np.abs(correction)) if not np.all(np.isfinite(correction)): raise RuntimeError('Nonfinite Newton correction; check initial guess.') if change <= np.finfo(float).eps*max(1.0, np.max(np.abs(y))): raise RuntimeError('Newton iteration cannot advance; check tolerance and initial guess.') y[1:n] += correction raise RuntimeError('Newton iteration did not converge within the iteration limit.') 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()