#! /usr/bin/env python3 # def humps_solve_bvp ( ): #*****************************************************************************80 # ## humps_solve_bvp() solves a boundary value problem. # # Discussion: # # The problem is posed on the interval [0,2]. # The exact solution is the humps function. # The software requires reformulating the second order BVP as a # pair of first order equations. # from humps import humps_fun from scipy.integrate import solve_bvp import matplotlib.pyplot as plt import numpy as np print ( '' ) print ( 'humps_bvp2():' ) print ( ' Solve a BVP using scipy.integrate.solve_bvp().' ) xa = 0.0 xb = 2.0 n = 21 x = np.linspace ( xa, xb, n ) m = 2 y = np.zeros ( [ m, n ] ) sol = solve_bvp ( humps_bvp_rhs, humps_bc, x, y ) x = sol.x y = sol.y # # Plot the computed data points, and the exact solution curve. # plt.plot ( x, y[0,:], 'ro' ) x2 = np.linspace ( xa, xb, 101 ) y2 = humps_fun ( x2 ) plt.plot ( x2, y2, 'b-' ) plt.grid ( True ) filename = 'humps_solve_bvp.png' plt.savefig ( filename ) print ( ' Graphics saved as "' + filename + '"' ) plt.show ( ) return def humps_bc ( ya, yb ): #*****************************************************************************80 # ## humps_bc() returns the boundary condition residual. # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 March 2025 # # Author: # # John Burkardt # # Input: # # real ya, yb: the current estimates for the solution at the # left and right endpoints. # # Output: # # real resid[2]: the residual in the boundary conditions. # from humps import humps_fun import numpy as np xa = 0.0 xb = 2.0 resid = np.array ( [ ya[0] - humps_fun(xa), yb[0] - humps_fun(xb) ] ) return resid def humps_bvp_rhs ( x, y ): #*****************************************************************************80 # ## humps_bvp_rhs() evaluates the right hand side of the humps bvp. # # Discussion: # # scipy.integrate.solve_bvp() requires us to transform a second order BVP # y""=f(x,y) # into a pair of first order equations: # y0' = y1 # y1' = f(x,y) # # y(x) = 1.0 / ( ( x - 0.3 )**2 + 0.01 ) \ # + 1.0 / ( ( x - 0.9 )**2 + 0.04 ) \ # - 6.0 # # f(x,y) = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 \ # - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2 # # Licensing: # # This code is distributed under the MIT license. # # Modified: # # 18 March 2025 # # Author: # # John Burkardt # # Input: # # real x(n): the arguments. # # Output: # # real ypp(2,n): the value of the right hand side functions. # import numpy as np u1 = - 2.0 * ( x - 0.3 ) v1 = ( ( x - 0.3 )**2 + 0.01 )**2 u2 = - 2.0 * ( x - 0.9 ) v2 = ( ( x - 0.9 )**2 + 0.04 )**2 u1p = - 2.0 v1p = 2.0 * ( ( x - 0.3 )**2 + 0.01 ) * 2.0 * ( x - 0.3 ) u2p = - 2.0 v2p = 2.0 * ( ( x - 0.9 )**2 + 0.04 ) * 2.0 * ( x - 0.9 ) n = len ( x ) rhs = np.zeros ( [ 2, n ] ) rhs[0] = y[1] rhs[1] = ( u1p * v1 - u1 * v1p ) / v1**2 \ + ( u2p * v2 - u2 * v2p ) / v2**2 return rhs if ( __name__ == "__main__" ): humps_solve_bvp ( )