|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 |
import os import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import sympy from IPython.display import display sympy.init_printing() mpl.rcParams['text.usetex'] = True import sympy os.chdir(r'D:\projects\wordpress\ex33') os.getcwd() def plot_direction_field(x, y_x, f_xy, x_lim=(-5, 5), y_lim=(-5, 5), ax=None): f_np = sympy.lambdify((x, y_x), f_xy, 'numpy') x_vec = np.linspace(x_lim[0], x_lim[1], 20) y_vec = np.linspace(y_lim[0], y_lim[1], 20) if ax is None: _, ax = plt.subplots(figsize=(4, 4)) dx = x_vec[1] - x_vec[0] dy = y_vec[1] - y_vec[0] for m, xx in enumerate(x_vec): for n, yy in enumerate(y_vec): Dy = f_np(xx, yy) * dx Dx = 0.8 * dx**2 / np.sqrt(dx**2 + Dy**2) Dy = 0.8 * Dy*dy / np.sqrt(dx**2 + Dy**2) ax.plot([xx - Dx/2, xx + Dx/2],[yy - Dy/2, yy + Dy/2], 'b', lw=0.5) ax.axis('tight') ax.set_title(r"$%s$" % (sympy.latex(sympy.Eq(y(x).diff(x), f_xy))), fontsize=18) return ax def apply_ics(sol, ics, x, known_params): free_params = sol.free_symbols - set(known_params) eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics) for n in range(len(ics))] sol_params = sympy.solve(eqs, free_params) return sol.subs(sol_params) x = sympy.symbols("x") y = sympy.Function("y") f = y(x)**2 + x display(sympy.Eq(y(x).diff(x), f)) ics = {y(0): 0} ode_sol = sympy.dsolve(y(x).diff(x) - f) display(ode_sol) ode_sol = apply_ics(ode_sol, {y(0): 0}, x, []) display(ode_sol) ode_sol = sympy.dsolve(y(x).diff(x) - f, ics=ics) display(ode_sol) fig, axes = plt.subplots(1, 2, figsize=(8, 4)) plot_direction_field(x, y(x), f, ax=axes[0]) x_vec = np.linspace(-3, 3, 100) axes[0].plot(x_vec, sympy.lambdify(x, ode_sol.rhs.removeO())(x_vec), 'b', lw=2) axes[0].set_ylim(-5, 5) plot_direction_field(x, y(x), f, ax=axes[1]) x_vec = np.linspace(-1, 1, 100) axes[1].plot(x_vec, sympy.lambdify(x, ode_sol.rhs.removeO())(x_vec), 'b', lw=2) ode_sol_m = ode_sol_p = ode_sol dx = 0.125 for x0 in np.arange(1, 2., dx): x_vec = np.linspace(x0, x0 + dx, 100) ics = {y(x0): ode_sol_p.rhs.removeO().subs(x, x0)} ode_sol_p = sympy.dsolve(y(x).diff(x) - f, ics=ics, n=6) axes[1].plot(x_vec, sympy.lambdify(x, ode_sol_p.rhs.removeO())(x_vec), 'r', lw=2) for x0 in np.arange(1, 5, dx): x_vec = np.linspace(-x0-dx, -x0, 100) ics = {y(-x0): ode_sol_m.rhs.removeO().subs(x, -x0)} ode_sol_m = sympy.dsolve(y(x).diff(x) - f, ics=ics, n=6) axes[1].plot(x_vec, sympy.lambdify(x, ode_sol_m.rhs.removeO())(x_vec), 'r', lw=2) fig.tight_layout() plt.savefig("example33.png", dpi=100) plt.show() plt.close() <img class=" wp-image-155 aligncenter" src="https://gantovnik.com/bio-tips/wp-content/uploads/2019/01/example33.png" alt="example33" width="614" height="307" /> |
Recent Comments