From 1edc2434b94d47c806567a15be8d531f462a1d3a Mon Sep 17 00:00:00 2001 From: Gerardo Marx Date: Mon, 29 Jun 2020 19:08:50 -0500 Subject: [PATCH] Numerical approximations with series --- readme.md | 8 ++++++++ scripts/program-02d.py | 12 ++++++++++++ scripts/program-02e.py | 29 +++++++++++++++++++++++++++++ 3 files changed, 49 insertions(+) create mode 100644 scripts/program-02d.py create mode 100644 scripts/program-02e.py diff --git a/readme.md b/readme.md index 5e19599..444a375 100644 --- a/readme.md +++ b/readme.md @@ -1,2 +1,10 @@ # Dynamical systems examples Lynch Repository for examples related with book Dynamical system with applications using Python by Stephen Lynch. + +# List of examples included +## Chapter 2 +- program 2-a: +- program 2-b: +- program 2-c: +- program 2-d: +- program 2-e: Numerical and truncated series solutions. diff --git a/scripts/program-02d.py b/scripts/program-02d.py new file mode 100644 index 0000000..f6c05ac --- /dev/null +++ b/scripts/program-02d.py @@ -0,0 +1,12 @@ +# Program 02d: Power series solution of a second order ODE; +# example 8, chapter 2 +from sympy import Function, dsolve, pprint +from sympy.abc import t + +x = Function('x') +dxdt2 = x(t).diff(t, 2) +dxdt = x(t).diff(t) +ode = dxdt2+2*t**2*dxdt+x(t) +sol = dsolve(ode, hint='2nd_power_series_ordinary', n=6) +print('Solution:\n') +pprint(sol) diff --git a/scripts/program-02e.py b/scripts/program-02e.py new file mode 100644 index 0000000..e3d843a --- /dev/null +++ b/scripts/program-02e.py @@ -0,0 +1,29 @@ +# Program 02e: Numerical and truncated series solution +# figure 2.6 +import matplotlib.pyplot as plt +import numpy as np +from scipy.integrate import odeint + + +def ode(X, t): + x = X[0] + y = X[1] + dxdt = y + dydt = x-t**2*y + return [dxdt, dydt] + + +X0 = [1, 0] +t = np.linspace(0, 10, 1000) +sol = odeint(ode, X0, t) +x = sol[:, 0] +y = sol[:, 1] +fig, ax = plt.subplots() +ax.plot(t, x, label='Numerical') +ax.plot(t, 1+t**2/2+t**4/24, 'r-', label='Truncated series') +plt.xlabel('t') +plt.ylabel('x') +plt.xlim(0, 4) +plt.ylim(0, 4) +ax.legend() +plt.show()