|
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 |
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() #%matplotlib inline #%config InlineBackend.figure_format='retina' os.chdir(r'D:\projects\wordpress\ex30') os.getcwd() #Symbolic ODE solving with SymPy #Newton's law of cooling t, k, T0, Ta = sympy.symbols("t, k, T_0, T_a") T = sympy.Function("T") ode = T(t).diff(t) + k*(T(t) - Ta) ode_sol = sympy.dsolve(ode) display(ode_sol.lhs) display(ode_sol.rhs) ics = {T(0): T0} display(ics) C_eq = sympy.Eq(ode_sol.lhs.subs(t, 0).subs(ics), ode_sol.rhs.subs(t, 0)) display(C_eq) C_sol = sympy.solve(C_eq) display(C_sol) display(ode_sol.subs(C_sol[0]))</pre> def apply_ics(sol, ics, x, known_params): # Apply the initial conditions (ics), given as a dictionary on # the form ics = {y(0): y0: y(x).diff(x).subs(x, 0): yp0, ...} # to the solution of the ODE with indepdendent variable x. # The undetermined integration constants C1, C2, ... are extracted # from the free symbols of the ODE solution, excluding symbols in # the known_params list. 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) display(ode_sol) apply_ics(ode_sol, ics, t, [k, Ta]) ode_sol = apply_ics(ode_sol, ics, t, [k, Ta]).simplify() display(ode_sol) y_x = sympy.lambdify((t, k), ode_sol.rhs.subs({T0: 5, Ta: 1}), 'numpy') fig, ax = plt.subplots(figsize=(8, 4)) x = np.linspace(0, 4, 100) for k in [1, 2, 3]: ax.plot(x, y_x(x, k), label=r"$k=%d$" % k) ax.set_title(r"$%s$" % sympy.latex(ode_sol), fontsize=18) ax.set_xlabel(r"$x$", fontsize=18) ax.set_ylabel(r"$y$", fontsize=18) ax.legend() fig.tight_layout() plt.savefig("example30.png", dpi=100) plt.show() plt.close() |
Recent Comments