• API
• FAQ
• Tools
• Archive
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, interval, N+2)
28.     y = np.zeros_like(s)
29.     y = bcs
30.
31.     def rhs_shooting(s, z):
32.         """
33.         RHS function for shooting (z = [y, y']; solve z' = [z, 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, f(x, z, z)]
48.         """
49.
50.         v = z
51.         theta = z
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
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, guess])
74.         return sol_shoot.y[1, -1] - bcs
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, 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.

Top