Демонстрация. Методы научного программирования на примере библиотеки символьных вычислений SymPy.¶

На основе урока от Jülich Supercomputing Centre.

Добавляйте ячейки исходного кода программы:

In [1]:

from IPython.display import display

from sympy.interactive import printing
printing.init_printing()

from __future__ import division
import sympy as sym
from sympy import *
x, y, z = symbols("x y z")
k, m, n = symbols("k m n", integer=True)
f, g, h = map(Function, 'fgh')

from IPython.display import display from sympy.interactive import printing printing.init_printing() from __future__ import division import sympy as sym from sympy import * x, y, z = symbols("x y z") k, m, n = symbols("k m n", integer=True) f, g, h = map(Function, 'fgh')

Добавляйте подзаголовки

И свои комментарии:

Исполняемые фрагменты кода:

Алгебра

In [7]:

eq = ((x+y)**2 * (x+1))
eq

eq = ((x+y)**2 * (x+1)) eq

Out[7]:

$$\left(x + 1\right) \left(x + y\right)^{2}$$

In [8]:

expand(eq)

expand(eq)

Out[8]:

$$x^{3} + 2 x^{2} y + x^{2} + x y^{2} + 2 x y + y^{2}$$

In [9]:

a = 1/x + (x*sin(x) - 1)/x
a

a = 1/x + (x*sin(x) - 1)/x a

Out[9]:

$$\frac{1}{x} \left(x \sin{\left (x \right )} - 1\right) + \frac{1}{x}$$

In [10]:

simplify(a)

simplify(a)

Out[10]:

$$\sin{\left (x \right )}$$

In [11]:

eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0)
eq

eq = Eq(x**3 + 2*x**2 + 4*x + 8, 0) eq

Out[11]:

$$x^{3} + 2 x^{2} + 4 x + 8 = 0$$

In [12]:

solve(eq, x)

solve(eq, x)

Out[12]:

$$\begin{bmatrix}-2, & - 2 i, & 2 i\end{bmatrix}$$

In [13]:

a, b = symbols('a b')
Sum(6*n**2 + 2**n, (n, a, b))

a, b = symbols('a b') Sum(6*n**2 + 2**n, (n, a, b))

Out[13]:

$$\sum_{n=a}^{b} \left(2^{n} + 6 n^{2}\right)$$

Мат. анализ

In [14]:

limit((sin(x)-x)/x**3, x, 0)

limit((sin(x)-x)/x**3, x, 0)

Out[14]:

$$- \frac{1}{6}$$

In [15]:

(1/cos(x)).series(x, 0, 6)

(1/cos(x)).series(x, 0, 6)

Out[15]:

$$1 + \frac{x^{2}}{2} + \frac{5 x^{4}}{24} + \mathcal{O}\left(x^{6}\right)$$

In [16]:

diff(cos(x**2)**2 / (1+x), x)

diff(cos(x**2)**2 / (1+x), x)

Out[16]:

$$- \frac{4 x \cos{\left (x^{2} \right )}}{x + 1} \sin{\left (x^{2} \right )} - \frac{\cos^{2}{\left (x^{2} \right )}}{\left(x + 1\right)^{2}}$$

In [17]:

integrate(x**2 * cos(x), (x, 0, pi/2))

integrate(x**2 * cos(x), (x, 0, pi/2))

Out[17]:

$$-2 + \frac{\pi^{2}}{4}$$

In [18]:

eqn = Eq(Derivative(f(x),x,x) + 9*f(x), 1)
display(eqn)
dsolve(eqn, f(x))

eqn = Eq(Derivative(f(x),x,x) + 9*f(x), 1) display(eqn) dsolve(eqn, f(x))

$$9 f{\left (x \right )} + \frac{d^{2}}{d x^{2}} f{\left (x \right )} = 1$$

Out[18]:

$$f{\left (x \right )} = C_{1} \sin{\left (3 x \right )} + C_{2} \cos{\left (3 x \right )} + \frac{1}{9}$$

Научная визуализация¶

На примере разложения в ряд тейлора.

In [19]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline import numpy as np import matplotlib.pyplot as plt

In [21]:

def plot_taylor_approximations(func, x0=None, orders=(2, 4), xrange=(0,1), yrange=None, npts=200):
    """Plot the Taylor series approximations to a function at various orders.

    Parameters
    ----------
    func : a sympy function
    x0 : float
      Origin of the Taylor series expansion.  If not given, x0=xrange[0].
    orders : list
      List of integers with the orders of Taylor series to show.  Default is (2, 4).
    xrange : 2-tuple or array.
      Either an (xmin, xmax) tuple indicating the x range for the plot (default is (0, 1)),
      or the actual array of values to use.
    yrange : 2-tuple
      (ymin, ymax) tuple indicating the y range for the plot.  If not given,
      the full range of values will be automatically used. 
    npts : int
      Number of points to sample the x range with.  Default is 200.
    """
    if not callable(func):
        raise ValueError('func must be callable')
    if isinstance(xrange, (list, tuple)):
        x = np.linspace(float(xrange[0]), float(xrange[1]), npts)
    else:
        x = xrange
    if x0 is None: x0 = x[0]
    xs = sym.Symbol('x')
    # Make a numpy-callable form of the original function for plotting
    fx = func(xs)
    f = sym.lambdify(xs, fx, modules=['numpy'])
    # We could use latex(fx) instead of str(), but matploblib gets confused
    # with some of the (valid) latex constructs sympy emits.  So we play it safe.
    plt.plot(x, f(x), label=str(fx), lw=2)
    # Build the Taylor approximations, plotting as we go
    apps = {}
    for order in orders:
        app = fx.series(xs, x0, n=order).removeO()
        apps[order] = app
        # Must be careful here: if the approximation is a constant, we can't
        # blindly use lambdify as it won't do the right thing.  In that case, 
        # evaluate the number as a float and fill the y array with that value.
        if isinstance(app, sym.numbers.Number):
            y = np.zeros_like(x)
            y.fill(app.evalf())
        else:
            fa = sym.lambdify(xs, app, modules=['numpy'])
            y = fa(x)
        tex = sym.latex(app).replace('$', '')
        plt.plot(x, y, label=r'$n=%s:\, %s$' % (order, tex) )
        
    # Plot refinements
    if yrange is not None:
        plt.ylim(*yrange)
    plt.grid()
    plt.legend(loc='best').get_frame().set_alpha(0.8)

def plot_taylor_approximations(func, x0=None, orders=(2, 4), xrange=(0,1), yrange=None, npts=200): """Plot the Taylor series approximations to a function at various orders. Parameters ---------- func : a sympy function x0 : float Origin of the Taylor series expansion. If not given, x0=xrange[0]. orders : list List of integers with the orders of Taylor series to show. Default is (2, 4). xrange : 2-tuple or array. Either an (xmin, xmax) tuple indicating the x range for the plot (default is (0, 1)), or the actual array of values to use. yrange : 2-tuple (ymin, ymax) tuple indicating the y range for the plot. If not given, the full range of values will be automatically used. npts : int Number of points to sample the x range with. Default is 200. """ if not callable(func): raise ValueError('func must be callable') if isinstance(xrange, (list, tuple)): x = np.linspace(float(xrange[0]), float(xrange[1]), npts) else: x = xrange if x0 is None: x0 = x[0] xs = sym.Symbol('x') # Make a numpy-callable form of the original function for plotting fx = func(xs) f = sym.lambdify(xs, fx, modules=['numpy']) # We could use latex(fx) instead of str(), but matploblib gets confused # with some of the (valid) latex constructs sympy emits. So we play it safe. plt.plot(x, f(x), label=str(fx), lw=2) # Build the Taylor approximations, plotting as we go apps = {} for order in orders: app = fx.series(xs, x0, n=order).removeO() apps[order] = app # Must be careful here: if the approximation is a constant, we can't # blindly use lambdify as it won't do the right thing. In that case, # evaluate the number as a float and fill the y array with that value. if isinstance(app, sym.numbers.Number): y = np.zeros_like(x) y.fill(app.evalf()) else: fa = sym.lambdify(xs, app, modules=['numpy']) y = fa(x) tex = sym.latex(app).replace('$', '') plt.plot(x, y, label=r'$n=%s:\, %s$' % (order, tex) ) # Plot refinements if yrange is not None: plt.ylim(*yrange) plt.grid() plt.legend(loc='best').get_frame().set_alpha(0.8)

Определив функцию, строим графики:

In [22]:

plot_taylor_approximations(sin, 0, [2, 4, 6], (0, 2*pi), (-2,2))

plot_taylor_approximations(sin, 0, [2, 4, 6], (0, 2*pi), (-2,2))

In [23]:

plot_taylor_approximations(cos, 0, [2, 4, 6], (0, 2*pi), (-2,2))

plot_taylor_approximations(cos, 0, [2, 4, 6], (0, 2*pi), (-2,2))

In [24]:

plot_taylor_approximations(lambda x: 1/cos(x), 0, [2,4,6], (0, 2*pi), (-5,5))

plot_taylor_approximations(lambda x: 1/cos(x), 0, [2,4,6], (0, 2*pi), (-5,5))

Комментарии