|
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 |
import os import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['text.usetex'] = True import sympy from IPython.display import display sympy.init_printing() 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) os.chdir(r'D:\projects\wordpress\ex31') os.getcwd() t, k, T0, Ta = sympy.symbols("t, k, T_0, T_a") t, omega0 = sympy.symbols("t, omega_0", positive=True) gamma = sympy.symbols("gamma", complex=True) x = sympy.Function("x") ode = x(t).diff(t, 2) + 2 * gamma * omega0 * x(t).diff(t) + omega0**2 * x(t) display(ode) ode_sol = sympy.dsolve(ode) display(ode_sol) ics = {x(0): 1, x(t).diff(t).subs(t, 0): 0} display(ics) x_t_sol = apply_ics(ode_sol, ics, t, [omega0, gamma]) display(x_t_sol) x_t_critical = sympy.limit(x_t_sol.rhs, gamma, 1) display(x_t_critical) fig, ax = plt.subplots(figsize=(8, 4)) tt = np.linspace(0, 3, 250) for g in [0.1, 0.5, 1, 2.0, 5.0]: if g == 1: x_t = sympy.lambdify(t, x_t_critical.subs({omega0: 2.0 * sympy.pi}), 'numpy') else: x_t = sympy.lambdify(t, x_t_sol.rhs.subs({omega0: 2.0 * sympy.pi, gamma: g}), 'numpy') ax.plot(tt, x_t(tt).real, label=r"$\gamma = %.1f$" % g) ax.set_xlabel(r"$t$", fontsize=18) ax.set_ylabel(r"$x(t)$", fontsize=18) ax.legend() fig.tight_layout() plt.savefig("example31.png", dpi=100) plt.show() plt.close() <img class=" wp-image-149 aligncenter" src="https://gantovnik.com/bio-tips/wp-content/uploads/2019/01/example31.png" alt="example31" width="588" height="294" /> |
Recent Comments