In [1]:
from sympy import *
init_printing(use_unicode=True)
import matplotlib.pyplot as plt
%matplotlib inline
Vector Calculus¶
In [2]:
from sympy.vector import *
In [3]:
#Define the coordinate system in 3D
R = CoordSys3D('R')
In [4]:
v = 3*R.i + 4*R.j + 5*R.k
v
Out[4]:
$$(3)\mathbf{\hat{i}_{R}} + (4)\mathbf{\hat{j}_{R}} + (5)\mathbf{\hat{k}_{R}}$$
In [5]:
v.components
Out[5]:
$$\left \{ \mathbf{\hat{i}_{R}} : 3, \quad \mathbf{\hat{j}_{R}} : 4, \quad \mathbf{\hat{k}_{R}} : 5\right \}$$
Cross¶
In [6]:
#Obtain perpendicular direction
R.i.cross(R.j)
Out[6]:
$$\mathbf{\hat{k}_{R}}$$
In [7]:
#or shorter form
R.i^R.i
Out[7]:
$$\mathbf{\hat{0}}$$
In [8]:
v^R.i
Out[8]:
$$(5)\mathbf{\hat{j}_{R}} + (-4)\mathbf{\hat{k}_{R}}$$
dot¶
In [9]:
R.i.dot(R.j)
Out[9]:
$$0$$
In [11]:
#Shorter form
R.i&R.i
Out[11]:
$$1$$
In [12]:
v.dot(R.k)
Out[12]:
$$5$$
In [13]:
v.magnitude()
Out[13]:
$$5 \sqrt{2}$$
In [14]:
v.normalize()
Out[14]:
$$(\frac{3 \sqrt{2}}{10})\mathbf{\hat{i}_{R}} + (\frac{2 \sqrt{2}}{5})\mathbf{\hat{j}_{R}} + (\frac{\sqrt{2}}{2})\mathbf{\hat{k}_{R}}$$
Curl¶
Rotacional
In [15]:
v1 = R.y*R.z*R.i + R.x*R.z*R.j + R.x*R.y*R.k
v1
Out[15]:
$$(\mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{i}_{R}} + (\mathbf{{x}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{j}_{R}} + (\mathbf{{x}_{R}} \mathbf{{y}_{R}})\mathbf{\hat{k}_{R}}$$
In [16]:
curl(v1)
Out[16]:
$$\mathbf{\hat{0}}$$
In [17]:
v2 = R.x*R.y*R.z*R.i
v2
Out[17]:
$$(\mathbf{{x}_{R}} \mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{i}_{R}}$$
In [18]:
curl(v2)
Out[18]:
$$(\mathbf{{x}_{R}} \mathbf{{y}_{R}})\mathbf{\hat{j}_{R}} + (- \mathbf{{x}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{k}_{R}}$$
In [19]:
v1 = R.i + R.j + R.k
curl(p*v1)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-19-46cf8ff55cdb> in <module>() 1 v1 = R.i + R.j + R.k ----> 2 curl(p*v1) NameError: name 'p' is not defined
Divergence¶
In [33]:
v1 = R.x*R.y*R.z * (R.i+R.j+R.k)
v1
Out[33]:
$$(\mathbf{{x}_{R}} \mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{i}_{R}} + (\mathbf{{x}_{R}} \mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{j}_{R}} + (\mathbf{{x}_{R}} \mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{k}_{R}}$$
In [32]:
divergence(v1)
Out[32]:
$$\mathbf{{x}_{R}} \mathbf{{y}_{R}} + \mathbf{{x}_{R}} \mathbf{{z}_{R}} + \mathbf{{y}_{R}} \mathbf{{z}_{R}}$$
In [34]:
v2 = 2*R.y*R.z*R.j
v2
Out[34]:
$$(2 \mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{j}_{R}}$$
In [35]:
divergence(v2)
Out[35]:
$$2 \mathbf{{z}_{R}}$$
GRadient¶
In [36]:
s1 = R.x*R.y*R.z
s1
Out[36]:
$$\mathbf{{x}_{R}} \mathbf{{y}_{R}} \mathbf{{z}_{R}}$$
In [37]:
gradient(s1)
Out[37]:
$$(\mathbf{{y}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{i}_{R}} + (\mathbf{{x}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{j}_{R}} + (\mathbf{{x}_{R}} \mathbf{{y}_{R}})\mathbf{\hat{k}_{R}}$$
Line integral¶
In [126]:
from sympy import Curve, line_integrate, E, ln
from sympy.abc import x, y, t
In [127]:
C = Curve([E**t + 1, E**t - 1], (t, 0, ln(2)))
C
Out[127]:
$$Curve((exp(t) + 1, exp(t) - 1), (t, 0, log(2)))$$
In [129]:
line_integrate(x +y, C, [x, y])
Out[129]:
$$3 \sqrt{2}$$
Spherical coord¶
x,y,z -> spherical
In [98]:
from sympy.vector import CoordSys3D
a = CoordSys3D('a')
In [102]:
#b = a.create_new('b', transformation='spherical')
In [103]:
#b.transformation_to_parent()
In [104]:
#b.transformation_from_parent()
Fourier¶
In [20]:
x, y, z = symbols('x y z')
In [21]:
s = fourier_series(x**2, (x, -pi, pi))
s
Out[21]:
$$- 4 \cos{\left (x \right )} + \cos{\left (2 x \right )} + \frac{\pi^{2}}{3} + \ldots$$
In [22]:
s.truncate(4)
Out[22]:
$$- 4 \cos{\left (x \right )} + \cos{\left (2 x \right )} - \frac{4}{9} \cos{\left (3 x \right )} + \frac{\pi^{2}}{3}$$
Fourier transform¶
In [112]:
k =symbols('k')
In [113]:
fourier_transform(exp(-x**2), x, k)
Out[113]:
$$\sqrt{\pi} e^{- \pi^{2} k^{2}}$$
In [114]:
inverse_fourier_transform(sqrt(pi)*exp(-(pi*k)**2), k, x)
Out[114]:
$$e^{- x^{2}}$$
Piecewise¶
In [23]:
f = x #(x+1)**2
g= 9 #log(x)**2
p = Piecewise((0, x <= -1), (f, x <= 3), (g, True))
In [24]:
p.subs(x, 1)
Out[24]:
$$1$$
In [25]:
i = [i for i in range(-10, 10)]
test = map(lambda i: p.subs(x, i).evalf(), i )
plt.plot(i, test)
plt.grid()
Heaviside¶
In [115]:
Heaviside(9)
Out[115]:
$$1$$
In [116]:
Heaviside(-9)
Out[116]:
$$0$$
In [117]:
Heaviside(0)
Out[117]:
$$\theta\left(0\right)$$
In [118]:
Heaviside(x).fdiff()
Out[118]:
$$\delta\left(x\right)$$
Laplace¶
In [124]:
from sympy.integrals import laplace_transform
from sympy.abc import t, s, a
In [125]:
laplace_transform(t**a, t, s)
Out[125]:
$$\left ( \frac{s^{- a}}{s} \Gamma{\left(a + 1 \right)}, \quad 0, \quad - \Re{\left(a\right)} < 1\right )$$
Diff eqn¶
In [133]:
f = Function('f')
In [134]:
eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x)
eq
Out[134]:
$$\sin{\left (x \right )} \cos{\left (f{\left (x \right )} \right )} + \sin{\left (f{\left (x \right )} \right )} \cos{\left (x \right )} \frac{d}{d x} f{\left (x \right )}$$
In [135]:
dsolve(eq, hint='1st_exact')
Out[135]:
$$\left [ f{\left (x \right )} = - \operatorname{acos}{\left (\frac{C_{1}}{\cos{\left (x \right )}} \right )} + 2 \pi, \quad f{\left (x \right )} = \operatorname{acos}{\left (\frac{C_{1}}{\cos{\left (x \right )}} \right )}\right ]$$
In [ ]:
In [147]:
eq = f(x).diff(x,x)-0*f(x).diff(x)-4*f(x) - x**2*sin(x)#exp(x)
eq
Out[147]:
$$- x^{2} \sin{\left (x \right )} - 4 f{\left (x \right )} + \frac{d^{2}}{d x^{2}} f{\left (x \right )}$$
In [148]:
dsolve(eq)
Out[148]:
$$f{\left (x \right )} = C_{1} e^{- 2 x} + C_{2} e^{2 x} - \frac{x^{2}}{5} \sin{\left (x \right )} - \frac{4 x}{25} \cos{\left (x \right )} - \frac{2}{125} \sin{\left (x \right )}$$
In [151]:
dsolve(f(x).diff(x, 4) + 2*f(x).diff(x, 3) -
2*f(x).diff(x, 2) - 6*f(x).diff(x) + 5*f(x), f(x),
hint='nth_linear_constant_coeff_homogeneous')
Out[151]:
$$f{\left (x \right )} = \left(C_{1} + C_{2} x\right) e^{x} + \left(C_{3} \sin{\left (x \right )} + C_{4} \cos{\left (x \right )}\right) e^{- 2 x}$$
Wronskian¶
In [165]:
from sympy.abc import rho, phi
f = Function('f')(x, y)
g1 = Function('g')(x, y)
g2 = x**2 + 3*y
f2 = x**3 + 3*y
In [166]:
hessian(f, (x, y), [g1, g2])
Out[166]:
$$\left[\begin{matrix}0 & 0 & \frac{\partial}{\partial x} g{\left (x,y \right )} & \frac{\partial}{\partial y} g{\left (x,y \right )}\\0 & 0 & 2 x & 3\\\frac{\partial}{\partial x} g{\left (x,y \right )} & 2 x & \frac{\partial^{2}}{\partial x^{2}} f{\left (x,y \right )} & \frac{\partial^{2}}{\partial x\partial y} f{\left (x,y \right )}\\\frac{\partial}{\partial y} g{\left (x,y \right )} & 3 & \frac{\partial^{2}}{\partial x\partial y} f{\left (x,y \right )} & \frac{\partial^{2}}{\partial y^{2}} f{\left (x,y \right )}\end{matrix}\right]$$
In [168]:
X = Matrix([rho*cos(phi), rho*sin(phi)])
Y = Matrix([rho, phi])
In [169]:
X.jacobian(Y)
Out[169]:
$$\left[\begin{matrix}\cos{\left (\phi \right )} & - \rho \sin{\left (\phi \right )}\\\sin{\left (\phi \right )} & \rho \cos{\left (\phi \right )}\end{matrix}\right]$$
In [175]:
f = Function('f')
g = Function('g')
Out[175]:
f
In [183]:
wronskian([x, x**2], x)
Out[183]:
$$x^{2}$$
In [185]:
from IPython.display import Image
Image("figures/SL.jpg")
Out[185]:
TEnsor functions¶
In [186]:
KroneckerDelta(1, 2)
Out[186]:
$$0$$
In [187]:
KroneckerDelta(3, 3)
Out[187]:
$$1$$
In [188]:
LeviCivita(1, 2, 3)
Out[188]:
$$1$$
In [189]:
LeviCivita(1, 3, 2)
Out[189]:
$$-1$$
In [190]:
LeviCivita(1, 2, 2)
Out[190]:
$$0$$
Binomial¶
$\left(\frac{n}{k}\right)= \frac{n!}{k!(n-k)!}$
In [49]:
n = Symbol('n', integer=True, positive=True)
In [50]:
binomial(15, 8)
Out[50]:
$$6435$$
In [51]:
for N in range(8):
print([binomial(N, i) for i in range(N + 1)])
[1] [1, 1] [1, 2, 1] [1, 3, 3, 1] [1, 4, 6, 4, 1] [1, 5, 10, 10, 5, 1] [1, 6, 15, 20, 15, 6, 1] [1, 7, 21, 35, 35, 21, 7, 1]
In [52]:
binomial(n, 4).expand(func=True)
Out[52]:
$$\frac{n^{4}}{24} - \frac{n^{3}}{4} + \frac{11 n^{2}}{24} - \frac{n}{4}$$
In [53]:
expand_func(binomial(n, 4))
Out[53]:
$$\frac{n}{24} \left(n - 3\right) \left(n - 2\right) \left(n - 1\right)$$
Factorial¶
(n - 2k)
In [54]:
factorial(7)
Out[54]:
$$5040$$
In [75]:
var('n')
Out[75]:
$$n$$
In [77]:
factorial2(n + 1)
Out[77]:
$$\left(n + 1\right)!!$$
In [62]:
factorial2(5), factorial2(6)
Out[62]:
$$\left ( 15, \quad 48\right )$$
Fibonancci¶
In [63]:
[fibonacci(x) for x in range(11)]
Out[63]:
$$\left [ 0, \quad 1, \quad 1, \quad 2, \quad 3, \quad 5, \quad 8, \quad 13, \quad 21, \quad 34, \quad 55\right ]$$
Tribonancci¶
instead of starting with two predetermined terms, the sequence starts with three predetermined terms and each term afterwards is the sum of the preceding three terms.
0, 0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149, 274, 504, 927, 1705, 3136, 5768, 10609, 19513, 35890, 66012,
In [69]:
[tribonacci(x) for x in range(11)]
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-69-8fd3055a0a90> in <module>() ----> 1 [tribonacci(x) for x in range(11)] NameError: name 'tribonacci' is not defined
In [71]:
[harmonic(n) for n in range(6)]
Out[71]:
$$\left [ 0, \quad 1, \quad \frac{3}{2}, \quad \frac{11}{6}, \quad \frac{25}{12}, \quad \frac{137}{60}\right ]$$
In [72]:
[harmonic(n, 2) for n in range(6)]
Out[72]:
$$\left [ 0, \quad 1, \quad \frac{5}{4}, \quad \frac{49}{36}, \quad \frac{205}{144}, \quad \frac{5269}{3600}\right ]$$
$\sum_{n=1}^\infty \frac{1}{n^2}$, or $\zeta(2)$
In [73]:
harmonic(oo, 2)
Out[73]:
$$\frac{\pi^{2}}{6}$$
In [78]:
expand_func(harmonic(n+4))
Out[78]:
$$\operatorname{harmonic}{\left (n \right )} + \frac{1}{n + 4} + \frac{1}{n + 3} + \frac{1}{n + 2} + \frac{1}{n + 1}$$
In [79]:
from IPython.display import Image
In [82]:
Image("figures/lucas.png")
Out[82]:
In [84]:
# The golden ratio
[lucas(x) for x in range(11)]
Out[84]:
$$\left [ 2, \quad 1, \quad 3, \quad 4, \quad 7, \quad 11, \quad 18, \quad 29, \quad 47, \quad 76, \quad 123\right ]$$
In [ ]: