スカラー場とベクトル場 CoordSys3D

CoordSys3Dを使う方法とReferenceFrameを使う方法があるらしいです。
前者はsympy 1.1.1とnablaで検索して出てきたdocumentを実行してみました。

In [1]:
from sympy import init_printing
init_printing()
from sympy.vector import CoordSys3D
In [2]:
R=CoordSys3D('R')
v=3*R.i+4*R.j+5*R.k
v
Out[2]:
$$(3)\mathbf{\hat{i}_{R}} + (4)\mathbf{\hat{j}_{R}} + (5)\mathbf{\hat{k}_{R}}$$

i,j,kでなく、x,y,zでもいいみたいですが、a,b,cはダメでした。

In [3]:
w=3*R.x+4*R.y+5*R.z
w
Out[3]:
$$3 \mathbf{{x}_{R}} + 4 \mathbf{{y}_{R}} + 5 \mathbf{{z}_{R}}$$
In [4]:
electric_potential=2*R.x**2*R.y
electric_potential
Out[4]:
$$2 \mathbf{{x}_{R}}^{2} \mathbf{{y}_{R}}$$
In [5]:
from sympy import diff
diff(electric_potential,R.x)
Out[5]:
$$4 \mathbf{{x}_{R}} \mathbf{{y}_{R}}$$
In [6]:
v=R.x**2*R.i+2*R.x*R.z*R.k
v
Out[6]:
$$(\mathbf{{x}_{R}}^{2})\mathbf{\hat{i}_{R}} + (2 \mathbf{{x}_{R}} \mathbf{{z}_{R}})\mathbf{\hat{k}_{R}}$$

$\nabla{}$ nabla記号をDel()で代用できます。doit()で実行します。

In [7]:
from sympy.vector import CoordSys3D, Del
C=CoordSys3D('C')
delop=Del()
gradient_field=delop(C.x*C.y*C.z)
gradient_field
Out[7]:
$$(\frac{\partial}{\partial \mathbf{{x}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{i}_{C}} + (\frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{j}_{C}} + (\frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{k}_{C}}$$
In [8]:
gradient_field.doit()
Out[8]:
$$(\mathbf{{y}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{i}_{C}} + (\mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{j}_{C}} + (\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{k}_{C}}$$

ベクトル微分　回転curl

$\mathbf{F}=\left(F_x,F_y,F_z\right)$として
$\nabla{}\times{}\mathbf{F}= \left(\frac{\partial{}F_z}{\partial{}y}-\frac{\partial{}F_z}{\partial{}y}\right)\hat{\mathbf{i}}+ \left(\frac{\partial{}F_x}{\partial{}z}-\frac{\partial{}F_z}{\partial{}x}\right)\hat{\mathbf{j}}+ \left(\frac{\partial{}F_y}{\partial{}x}-\frac{\partial{}F_x}{\partial{}y}\right)\hat{\mathbf{k}} $

In [9]:
from sympy.vector import CoordSys3D,Del
C=CoordSys3D('C')
delop=Del()

delopを$\times{}$ベクトル外積で計算する方法
delopを^排他的論理和XORで計算する方法
delopをcurlで計算する方法
の３つがあります。
$\nabla{}\times{}(xyz,0,0)$の計算をしてみます

In [10]:
delop.cross(C.x*C.y*C.z*C.i)
Out[10]:
$$(\frac{d}{d \mathbf{{y}_{C}}} 0 - \frac{d}{d \mathbf{{z}_{C}}} 0)\mathbf{\hat{i}_{C}} + (- \frac{d}{d \mathbf{{x}_{C}}} 0 + \frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{j}_{C}} + (\frac{d}{d \mathbf{{x}_{C}}} 0 - \frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{k}_{C}}$$
In [11]:
delop.cross(C.x*C.y*C.z*C.i).doit()
Out[11]:
$$(\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{j}_{C}} + (- \mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{k}_{C}}$$
In [12]:
delop^C.x*C.y*C.z*C.i
Out[12]:
$$(\frac{d}{d \mathbf{{y}_{C}}} 0 - \frac{d}{d \mathbf{{z}_{C}}} 0)\mathbf{\hat{i}_{C}} + (- \frac{d}{d \mathbf{{x}_{C}}} 0 + \frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{j}_{C}} + (\frac{d}{d \mathbf{{x}_{C}}} 0 - \frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{k}_{C}}$$
In [13]:
(delop^C.x*C.y*C.z*C.i).doit()
Out[13]:
$$(\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{j}_{C}} + (- \mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{k}_{C}}$$
In [14]:
from sympy.vector import curl
curl(C.x*C.y*C.z*C.i)
Out[14]:
$$(\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{j}_{C}} + (- \mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{k}_{C}}$$

ベクトル微分 発散divergence

$\mathbf{F}=\left(U,V,W\right)$として
$\nabla{}\dot{}\mathbf{F}= \frac{\partial{}U}{\partial{}x}+\frac{\partial{}V}{\partial{}y}+\frac{\partial{}W}{\partial{}z} $

In [15]:
from sympy.vector import CoordSys3D,Del
C=CoordSys3D('C')
delop=Del()

$\mathbf{F}=\left(xyz,xyz,xyz\right)$として
$\nabla{}\dot{}\mathbf{F}=\nabla{}\dot{}\left(xyz,xyz,xyz\right)$の計算をしてみます
delopをベクトルの内積で計算する方法
delopを理論積&で計算する方法
delopをdivergenceで計算する方法
の3つの方法があります。

In [16]:
delop.dot(C.x*C.y*C.z*(C.i+C.j+C.k))
Out[16]:
$$\frac{\partial}{\partial \mathbf{{x}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right) + \frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right) + \frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right)$$
In [17]:
delop.dot(C.x*C.y*C.z*(C.i+C.j+C.k)).doit()
Out[17]:
$$\mathbf{{x}_{C}} \mathbf{{y}_{C}} + \mathbf{{x}_{C}} \mathbf{{z}_{C}} + \mathbf{{y}_{C}} \mathbf{{z}_{C}}$$
In [18]:
(delop & C.x*C.y*C.z*(C.i+C.j+C.k))
Out[18]:
$$\frac{\partial}{\partial \mathbf{{x}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right) + \frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right) + \frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right)$$
In [19]:
(delop & C.x*C.y*C.z*(C.i+C.j+C.k)).doit()
Out[19]:
$$\mathbf{{x}_{C}} \mathbf{{y}_{C}} + \mathbf{{x}_{C}} \mathbf{{z}_{C}} + \mathbf{{y}_{C}} \mathbf{{z}_{C}}$$
In [20]:
from sympy.vector import divergence
divergence(C.x*C.y*C.z*(C.i+C.j+C.k))
Out[20]:
$$\mathbf{{x}_{C}} \mathbf{{y}_{C}} + \mathbf{{x}_{C}} \mathbf{{z}_{C}} + \mathbf{{y}_{C}} \mathbf{{z}_{C}}$$

ベクトル微分　勾配gradient

$F=f\left(x,y,z\right)$として
$\nabla{}f=\frac{\partial{}f}{\partial{}x}\mathbf{\hat{i}}+ \frac{\partial{}f}{\partial{}y}\mathbf{\hat{j}}+ \frac{\partial{}f}{\partial{}z}\mathbf{\hat{k}} $

In [21]:
from sympy.vector import CoordSys3D,Del
C=CoordSys3D('C')
delop=Del()

計算方法は３つあります。
$f=xyz$で計算してみます。

In [22]:
delop.gradient(C.x*C.y*C.z)
Out[22]:
$$(\frac{\partial}{\partial \mathbf{{x}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{i}_{C}} + (\frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{j}_{C}} + (\frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{k}_{C}}$$
In [23]:
delop.gradient(C.x*C.y*C.z).doit()
Out[23]:
$$(\mathbf{{y}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{i}_{C}} + (\mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{j}_{C}} + (\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{k}_{C}}$$
In [24]:
delop(C.x*C.y*C.z)
Out[24]:
$$(\frac{\partial}{\partial \mathbf{{x}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{i}_{C}} + (\frac{\partial}{\partial \mathbf{{y}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{j}_{C}} + (\frac{\partial}{\partial \mathbf{{z}_{C}}}\left(\mathbf{{x}_{C}} \mathbf{{y}_{C}} \mathbf{{z}_{C}}\right))\mathbf{\hat{k}_{C}}$$
In [25]:
delop(C.x*C.y*C.z).doit()
Out[25]:
$$(\mathbf{{y}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{i}_{C}} + (\mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{j}_{C}} + (\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{k}_{C}}$$
In [26]:
from sympy.vector import gradient
gradient(C.x*C.y*C.z)
Out[26]:
$$(\mathbf{{y}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{i}_{C}} + (\mathbf{{x}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{j}_{C}} + (\mathbf{{x}_{C}} \mathbf{{y}_{C}})\mathbf{\hat{k}_{C}}$$

方向微分

$\left(\vec{\upsilon}\dot{}\nabla{}\right)\mathbf{F}\left(x\right)$

In [27]:
from sympy.vector import CoordSys3D,Del
C=CoordSys3D('C')
delop=Del()
vel=C.i+C.j+C.k
scalar_field=C.x*C.y*C.z
vector_field=C.x*C.y*C.z*C.i
In [28]:
(vel.dot(delop))(scalar_field)
Out[28]:
$$\mathbf{{x}_{C}} \mathbf{{y}_{C}} + \mathbf{{x}_{C}} \mathbf{{z}_{C}} + \mathbf{{y}_{C}} \mathbf{{z}_{C}}$$
In [29]:
(vel & delop)(vector_field)
Out[29]:
$$(\mathbf{{x}_{C}} \mathbf{{y}_{C}} + \mathbf{{x}_{C}} \mathbf{{z}_{C}} + \mathbf{{y}_{C}} \mathbf{{z}_{C}})\mathbf{\hat{i}_{C}}$$
In [30]:
from sympy.vector import directional_derivative
directional_derivative(C.x*C.y*C.z,3*C.i+4*C.j+C.k)
Out[30]:
$$\mathbf{{x}_{C}} \mathbf{{y}_{C}} + 4 \mathbf{{x}_{C}} \mathbf{{z}_{C}} + 3 \mathbf{{y}_{C}} \mathbf{{z}_{C}}$$

保存場とソレノイド場

$\mathbf{F}=\left(yz,xz,xy\right)$でいろいろ計算してみます。

In [31]:
from sympy.vector import CoordSys3D, is_conservative
R=CoordSys3D('R')
field=R.y*R.z*R.i+R.x*R.z*R.j+R.x*R.y*R.k
field
Out[31]:
$$(\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 [32]:
is_conservative(field)
Out[32]:
True
In [33]:
from sympy.vector import curl
curl(field)
Out[33]:
$$\mathbf{\hat{0}}$$
In [34]:
from sympy.vector import CoordSys3D, is_solenoidal
R=CoordSys3D('R')
field=R.y*R.z*R.i+R.x*R.z*R.j+R.x*R.y*R.k
field
Out[34]:
$$(\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 [35]:
is_solenoidal(field)
Out[35]:
True
In [36]:
from sympy.vector import divergence
divergence(field)
Out[36]:
$$0$$

スカラーポテンシャル

In [37]:
from sympy.vector import CoordSys3D, Point
from sympy.vector import scalar_potential_difference
R=CoordSys3D('R')
P=R.origin.locate_new('P',1*R.i+2*R.j+3*R.k)
P
Out[37]:
$$P$$
In [38]:
vectorfield=4*R.x*R.y*R.i+2*R.x**2*R.j
vectorfield
Out[38]:
$$(4 \mathbf{{x}_{R}} \mathbf{{y}_{R}})\mathbf{\hat{i}_{R}} + (2 \mathbf{{x}_{R}}^{2})\mathbf{\hat{j}_{R}}$$
In [39]:
scalar_potential_difference(vectorfield,R,R.origin, P)
Out[39]:
$$4$$

スカラー場とベクトル場 ReferenceFrame

基本的にCoordSys3Dと同じですので、章立てせず、ササっと流していきます。
CoordSys3DだとR.x、R.y、R.zですが、ReferenceFrameだとR[0]、R[1]、R[2]となるようです。

In [40]:
from sympy import init_printing
init_printing()
from sympy.physics.vector import ReferenceFrame
R=ReferenceFrame('R')
electric_potential=2*R[0]**2*R[1]
electric_potential
Out[40]:
$$2 R_{x}^{2} R_{y}$$

$F=x^2y$のとき
$\frac{\partial{}F}{\partial{}x}=4xy$
を解きます。

In [41]:
from sympy import diff
diff(electric_potential,R[0])
Out[41]:
$$4 R_{x} R_{y}$$
In [42]:
from sympy.physics.vector import dynamicsymbols,express
q=dynamicsymbols('q')
q
Out[42]:
$$q{\left (t \right )}$$

たぶんz軸方向に$q(t)$回転させた空間を設計し、electric_potentialを表現

In [43]:
R1=R.orientnew('R1',rot_type='Axis',amounts=[q,R.z])
In [44]:
express(electric_potential,R1,variables=True)
Out[44]:
$$2 \left(R_{1 x} \sin{\left (q{\left (t \right )} \right )} + R_{1 y} \cos{\left (q{\left (t \right )} \right )}\right) \left(R_{1 x} \cos{\left (q{\left (t \right )} \right )} - R_{1 y} \sin{\left (q{\left (t \right )} \right )}\right)^{2}$$

$R_x$が$R_{1x}\sin{}(q(t))+R_{1y}\cos{}(q(t))$に
$R_y$が$R_{1x}\cos{}(q(t))-R_{1y}\sin{}(q(t))$になり
$R_z$は$R_{1z}$のままとなっています。

In [45]:
from sympy.physics.vector import time_derivative
time_derivative(electric_potential,R)
Out[45]:
$$0$$
In [46]:
time_derivative(electric_potential,R1).simplify()
Out[46]:
$$2 \left(R_{1 x} \cos{\left (q{\left (t \right )} \right )} - R_{1 y} \sin{\left (q{\left (t \right )} \right )}\right) \left(\frac{3}{2} R_{1 x}^{2} \cos{\left (2 q{\left (t \right )} \right )} - \frac{R_{1 x}^{2}}{2} - 3 R_{1 x} R_{1 y} \sin{\left (2 q{\left (t \right )} \right )} - \frac{3}{2} R_{1 y}^{2} \cos{\left (2 q{\left (t \right )} \right )} - \frac{R_{1 y}^{2}}{2}\right) \frac{d}{d t} q{\left (t \right )}$$

回転curl 発散divergence 勾配gradient

CoordSys3Dのcurl、divergence、gradientと同じですが、 ReferenceFrameのRをつける必要があります。

In [47]:
from sympy.physics.vector import ReferenceFrame
R=ReferenceFrame('R')
from sympy.physics.vector import curl

$\mathbf{F}=\left(xyz,0,0\right)$として
$\nabla{}\times{}\mathbf{F}= \left(0,xy,-xz\right)$

In [48]:
field=R[0]*R[1]*R[2]*R.x
field
Out[48]:
$$R_{x} R_{y} R_{z}\mathbf{\hat{r}_x}$$
In [49]:
curl(field,R) # Rをつけないといけない
Out[49]:
$$R_{x} R_{y}\mathbf{\hat{r}_y} - R_{x} R_{z}\mathbf{\hat{r}_z}$$

$\mathbf{F}=\left(xyz,xyz,xyz\right)$として
$\nabla{}\dot{}\mathbf{F}=xy+yz+zx$

In [50]:
from sympy.physics.vector import divergence
field=R[0]*R[1]*R[2]*(R.x+R.y+R.z)
field
Out[50]:
$$R_{x} R_{y} R_{z}\mathbf{\hat{r}_x} + R_{x} R_{y} R_{z}\mathbf{\hat{r}_y} + R_{x} R_{y} R_{z}\mathbf{\hat{r}_z}$$
In [51]:
divergence(field,R)
Out[51]:
$$R_{x} R_{y} + R_{x} R_{z} + R_{y} R_{z}$$

$F=xyz$として
$\nabla{}F=\left(yz,zx,xz\right)$

In [52]:
from sympy.physics.vector import gradient
scalar_field=R[0]*R[1]*R[2]
scalar_field
Out[52]:
$$R_{x} R_{y} R_{z}$$
In [53]:
gradient(scalar_field,R)
Out[53]:
$$R_{y} R_{z}\mathbf{\hat{r}_x} + R_{x} R_{z}\mathbf{\hat{r}_y} + R_{x} R_{y}\mathbf{\hat{r}_z}$$

保存場 ソレノイド場 スカラーポテンシャル

In [54]:
from sympy.physics.vector import ReferenceFrame,is_conservative
R=ReferenceFrame('R')
field=R[1]*R[2]*R.x+R[0]*R[2]*R.y+R[0]*R[1]*R.z
field
Out[54]:
$$R_{y} R_{z}\mathbf{\hat{r}_x} + R_{x} R_{z}\mathbf{\hat{r}_y} + R_{x} R_{y}\mathbf{\hat{r}_z}$$
In [55]:
is_conservative(field)
Out[55]:
True
In [56]:
curl(field,R)
Out[56]:
$$0$$
In [57]:
from sympy.physics.vector import is_solenoidal
In [58]:
is_solenoidal(field)
Out[58]:
True
In [59]:
divergence(field,R)
Out[59]:
$$0$$
In [60]:
from sympy.physics.vector import scalar_potential
conservative_field=4*R[0]*R[1]*R[2]*R.x+2*R[0]**2*R[2]*R.y+2*R[0]**2*R[1]*R.z
conservative_field
Out[60]:
$$4 R_{x} R_{y} R_{z}\mathbf{\hat{r}_x} + 2 R_{x}^{2} R_{z}\mathbf{\hat{r}_y} + 2 R_{x}^{2} R_{y}\mathbf{\hat{r}_z}$$
In [61]:
scalar_potential(conservative_field,R)
Out[61]:
$$2 R_{x}^{2} R_{y} R_{z}$$
In [62]:
from sympy.physics.vector import Point,scalar_potential_difference
O=Point('O')
P=O.locatenew('P',1*R.x+2*R.y+3*R.z)
vector_field=4*R[0]*R[1]*R.x+2*R[0]**2*R.y
vector_field
Out[62]:
$$4 R_{x} R_{y}\mathbf{\hat{r}_x} + 2 R_{x}^{2}\mathbf{\hat{r}_y}$$
In [63]:
scalar_potential_difference(vector_field,R,O,P,O)
Out[63]:
$$4$$