SHARE
TWEET

Untitled

a guest Dec 5th, 2019 83 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. def shooting_dirichlet(f, interval, bcs, N):
  2.        
  3.     """
  4.     Solve a BVP using shooting
  5.    
  6.     Parameters
  7.     ----------
  8.    
  9.     f: function
  10.         RHS function defining BVP: y'' = f(x, y, y')
  11.     interval: list of two floats
  12.         [a, b] where BVP defined on x \in [a, b]
  13.     bcs: list of two floats
  14.         Defines Dirichlet boundary values [y(a), y(b)]
  15.     N: single integer
  16.         Number of interior points in the grid, used for output only
  17.    
  18.     Returns
  19.     -------
  20.    
  21.     x: array of floats
  22.         Location of grid points
  23.     y: array of floats
  24.         Solution at grid points
  25.     """
  26.    
  27.     s = np.linspace(interval[0], interval[1], N+2)
  28.     y = np.zeros_like(s)
  29.     y[0] = bcs[0]
  30.    
  31.     def rhs_shooting(s, z):
  32.         """
  33.         RHS function for shooting (z = [y, y']; solve z' = [z[1], f(...)])
  34.      
  35.         Parameters
  36.         ----------
  37.        
  38.         x: float
  39.         Location
  40.         z: list of floats
  41.         [y, y']
  42.        
  43.         Returns
  44.         -------
  45.        
  46.         dzdx: list of floats
  47.         [y', y''] = [z[1], f(x, z[0], z[1])]
  48.         """
  49.                
  50.         v = z[1]
  51.         theta = z[0]
  52.         dtheta = v
  53.         dvbyds = s * fg * np.cos(theta) + s * fx * np.sin(theta)
  54.         #y, dy = z
  55.         return np.array([dtheta, dvbyds])
  56.    
  57.     def residuals_shooting(guess):
  58.         """
  59.         Solve IVP given y'(a): compute r = y(b) - bcs[1]
  60.      
  61.         Parameters
  62.         ----------
  63.        
  64.         guess: float
  65.         Guess for derivative of solution at x=a
  66.      
  67.         Returns
  68.         -------
  69.        
  70.         residual: float
  71.         Values of residual r at x=b
  72.         """
  73.         sol_shoot = solve_ivp(rhs_shooting, interval, [bcs[0], guess])
  74.         return sol_shoot.y[1, -1] - bcs[1]
  75.    
  76.    
  77.     sol = scipy.optimize.root_scalar(residuals_shooting, method='brentq', bracket=[-10, 10])
  78.     guess = sol.root
  79.     sol_shoot = solve_ivp(rhs_shooting, interval, [bcs[0], guess], t_eval=s)
  80.     y_all = sol_shoot.y
  81.     return s, y_all[0, :]
  82.    
  83.         #Solve IVP with theta
  84. dxds = np.cos(y_all)
  85. dzds = np.sin(y_all)
  86. x0 = R*np.cos(theta0)
  87. z0 = R*np.sin(theta0)
  88.    
  89. x = solve_ivp(dxds, interval, x0)
  90. z = solve_ivp(dzds, interval, z0)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top