Return to SymPy で理工系の数学C

SymPy で偏微分

必要なパッケージを import します。

In [1]:
from sympy.abc import * 
from sympy import *

# SymPy Plotting Backends (SPB): グラフを描く際に利用
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']

偏微分：多変数関数の微分

SymPy での偏微分は（常微分と同様の書き方ですが）以下のように書きます。

$\displaystyle \frac{\partial}{\partial x} f(x, y) = $ diff(f(x, y), x);

$\displaystyle \frac{\partial}{\partial y} f(x, y) = $ diff(f(x, y), y);

In [2]:
f = Function('f')
In [3]:
diff(f(x, y), x)
Out[3]:
$\displaystyle \frac{\partial}{\partial x} f{\left(x,y \right)}$
In [4]:
diff(f(x, y), y)
Out[4]:
$\displaystyle \frac{\partial}{\partial y} f{\left(x,y \right)}$

例題

次の関数 $z$ の偏導関数を求めよ。

(1) $ z = x^3 – 4 x^2 y + x y + 3 y^2$

In [5]:
z = x**3 - 4*x**2*y + x*y + 3*y**2
In [6]:
diff(z, x)
Out[6]:
$\displaystyle 3 x^{2} – 8 x y + y$
In [7]:
diff(z, y)
Out[7]:
$\displaystyle – 4 x^{2} + x + 6 y$

(2) $\displaystyle z = \tan^{-1} \frac{y}{x} $

In [8]:
z = atan(y/x)
In [9]:
diff(z, x).simplify()
Out[9]:
$\displaystyle – \frac{y}{x^{2} + y^{2}}$
In [10]:
diff(z, y).simplify()
Out[10]:
$\displaystyle \frac{x}{x^{2} + y^{2}}$

全微分

2変数関数 $z = f(x, y)$ の全微分 $dz$ とは

$$ dz = df(x,y) = \frac{\partial f}{\partial x} dx + \frac{\partial f}{\partial y} dy$$

例題

次の関数 $z$ の全微分を求めよ。

(1) $ z = \sqrt{x^2 – y^2} $

(2) $ z = x y^2 $

In [11]:
z = sqrt(x**2 - y**2)
var('dx dy')

diff(z, x) * dx + diff(z, y) * dy
Out[11]:
$\displaystyle \frac{dx x}{\sqrt{x^{2} – y^{2}}} – \frac{dy y}{\sqrt{x^{2} – y^{2}}}$

上記の答え，
$\displaystyle dz = \frac{x}{\sqrt{x^{2} – y^{2}}} dx – \frac{y}{\sqrt{x^{2} – y^{2}}}dy $ と表示してほしいところ。

In [12]:
z = x*y**2
diff(z, x) * dx + diff(z, y) * dy
Out[12]:
$\displaystyle dx y^{2} + 2 dy x y$

上記の答え，
$\displaystyle dz = y^{2}\,dx + 2 x y\, dy $ と表示してほしいところ。

高階偏導関数

SymPy が2階偏導関数に関して

$$ \frac{\partial^2 f}{\partial x \partial y} = \frac{\partial^2 f}{\partial y \partial x} $$

のように，偏微分の順序を変えてもよいことを知っているか確認。

In [13]:
f = Function('f')

Eq((Derivative(Derivative(f(x, y), x), y) - 
    Derivative(Derivative(f(x, y), y), x)), 
   (Derivative(Derivative(f(x, y), x), y) - 
    Derivative(Derivative(f(x, y), y), x)).doit())
Out[13]:
$\displaystyle \frac{\partial^{2}}{\partial y\partial x} f{\left(x,y \right)} – \frac{\partial^{2}}{\partial x\partial y} f{\left(x,y \right)} = 0$

テイラー展開（2変数）

1変数関数 $f(x)$ のテイラー展開の復習

$f(x)$ を $x = a$ のまわりで $n = 5$ 次（$x^5$）まで展開するときは，SymPy では
series(f(x), x, a, 6) のように書く。

例として，$f(x) = e^{x}$ を $x = 0$ まわりで $5$ 次まで展開してみる。

In [14]:
f = exp(x)

series(f, x, 0, 6) # 5 次までなら 6 と書く。それが Python 流。
Out[14]:
$\displaystyle 1 + x + \frac{x^{2}}{2} + \frac{x^{3}}{6} + \frac{x^{4}}{24} + \frac{x^{5}}{120} + O\left(x^{6}\right)$

$\sin x$ を $x = 0$ のまわりで $x^3$ まで展開するときは…

In [15]:
series(sin(x), x, 0, 4) # 3 次までなら 4 と書く。
Out[15]:
$\displaystyle x – \frac{x^{3}}{6} + O\left(x^{4}\right)$

2変数関数 $f(x, y)$ のテイラー展開

本題である2変数関数 $f(x, y)$ のテイラー展開。series() は1変数関数のテイラー展開のみのようだ。

例として，$\displaystyle f(x, y) = \frac{e^x}{1-x+2 y}$ を $x = 0, y = 0$ のまわりで $2$ 次まで展開してみる。

In [16]:
f = exp(x)/(1-x+2*y)
f
Out[16]:
$\displaystyle \frac{e^{x}}{- x + 2 y + 1}$

$x \rightarrow \epsilon x, \ y \rightarrow \epsilon y$ とおき…

In [17]:
var('epsilon')
F = f.subs(x, epsilon*x).subs(y, epsilon*y)
F
Out[17]:
$\displaystyle \frac{e^{\epsilon x}}{- \epsilon x + 2 \epsilon y + 1}$

1.$\epsilon$ について2次までテイラー展開し，$\dots$ series()

2.カッコをはずして展開。$\dots$ expand()

In [18]:
series(F, epsilon, 0, 3).expand()
Out[18]:
$\displaystyle 1 – 2 \epsilon y + 2 \epsilon x + 4 \epsilon^{2} y^{2} – 6 \epsilon^{2} x y + \frac{5 \epsilon^{2} x^{2}}{2} + O\left(\epsilon^{3}\right)$

3.高次の項を示す $O()$ を取り去り，$\dots$ removeO()

（ここで降べきの順に並ぶ。）

In [19]:
_.removeO()
Out[19]:
$\displaystyle \frac{5 \epsilon^{2} x^{2}}{2} – 6 \epsilon^{2} x y + 4 \epsilon^{2} y^{2} + 2 \epsilon x – 2 \epsilon y + 1$

4.最後に $\epsilon \rightarrow 1$ とおきなおす。$\dots$ subs()

In [20]:
_.subs(epsilon, 1)
Out[20]:
$\displaystyle \frac{5 x^{2}}{2} – 6 x y + 2 x + 4 y^{2} – 2 y + 1$

合成関数の偏微分

例題

\begin{eqnarray}
z &=& f(x, y) = x y \\
x &=& x(u, v) = u\cos v \\
y &=& y(u, v) = u \sin v \\
\therefore\ \ z &=& f\left(x(u, v), y(u, v)\right) = z(u, v)
\end{eqnarray}

について，以下を求める。

$$ 1. \ \frac{\partial z}{\partial u} \quad 2. \ \frac{\partial z}{\partial v} $$

In [21]:
var('x y z u v')

x = u*cos(v)
y = u*sin(v)

z = x * y

$\displaystyle \frac{\partial z}{\partial u} = \cdots$

In [22]:
Eq(Derivative(z, u), 
   Derivative(z, u).doit())
Out[22]:
$\displaystyle \frac{\partial}{\partial u} u^{2} \sin{\left(v \right)} \cos{\left(v \right)} = 2 u \sin{\left(v \right)} \cos{\left(v \right)}$

右辺を簡単化してみます。

In [23]:
Eq(Derivative(z, u).doit(),
   Derivative(z, u).doit().simplify())
Out[23]:
$\displaystyle 2 u \sin{\left(v \right)} \cos{\left(v \right)} = u \sin{\left(2 v \right)}$

$\displaystyle \frac{\partial z}{\partial v} = \cdots$

In [24]:
Eq(Derivative(z, v), 
   Derivative(z, v).doit())
Out[24]:
$\displaystyle \frac{\partial}{\partial v} u^{2} \sin{\left(v \right)} \cos{\left(v \right)} = – u^{2} \sin^{2}{\left(v \right)} + u^{2} \cos^{2}{\left(v \right)}$
In [25]:
# 右辺を簡単化

Eq(diff(z, v), 
   diff(z, v).simplify())
Out[25]:
$\displaystyle – u^{2} \sin^{2}{\left(v \right)} + u^{2} \cos^{2}{\left(v \right)} = u^{2} \cos{\left(2 v \right)}$

陰関数定理

例題

$f(x, y) = x^2 + y^2 -1 = 0$ を満たす陰関数 $y=y(x)$ の微分を求める。

陰関数定理を使わずに解く

陰関数定理より，

$$ \frac{dy}{dx} = – \frac{\frac{\partial f}{\partial x}}{\frac{\partial f}{\partial y}}$$

まずはこれを使わずに SymPy で導いてみます。

まず，$f(x, y) = 0$ より $y$ は $x$ の（陰）関数となるので，$y(x)$。

In [26]:
var('x y f')

y = Function('y')
f = x**2 + y(x)**2 - 1
f
Out[26]:
$\displaystyle x^{2} + y^{2}{\left(x \right)} – 1$
In [27]:
df = diff(f, x)
df
Out[27]:
$\displaystyle 2 x + 2 y{\left(x \right)} \frac{d}{d x} y{\left(x \right)}$

$\displaystyle \frac{df}{dx} = 0$ を $\displaystyle \frac{dy}{dx}$ について解きます。

In [28]:
sol = solve(df, diff(y(x), x))

Eq(diff(y(x), x), sol[0])
Out[28]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = – \frac{x}{y{\left(x \right)}}$
陰関数定理を使って…

陰関数定理より，

$$ \frac{dy}{dx} = – \frac{\frac{\partial f}{\partial x}}{\frac{\partial f}{\partial y}}$$

In [29]:
var('x y f')
f = x**2 + y**2 - 1

dydx = - diff(f, x)/diff(f, y)
Eq(Derivative(y, x), dydx)
Out[29]:
$\displaystyle \frac{d}{d x} y = – \frac{x}{y}$
陽関数として表して微分する

別解として，陽関数として表して解く場合。

$f(x, y) = x^2 + y^2 – 1 = 0$ を $y$ について解く。

In [30]:
var('x y f')
f = x**2 + y**2 - 1
f
Out[30]:
$\displaystyle x^{2} + y^{2} – 1$
In [31]:
sols = solve(f, y)
sols
Out[31]:
[-sqrt(1 - x**2), sqrt(1 - x**2)]

1つ目の解は以下のように取り出す。

In [32]:
Eq(y, sols[0])
Out[32]:
$\displaystyle y = – \sqrt{1 – x^{2}}$

そしてそれを微分すると…

In [33]:
dy0 = diff(sols[0], x)
Eq(Derivative(y, x), dy0)
Out[33]:
$\displaystyle \frac{d}{d x} y = \frac{x}{\sqrt{1 – x^{2}}}$

微分した結果の右辺の分母に $\sqrt{1-x^2} = -y$ を代入して $y$ で書き直すと…

In [34]:
Eq(Derivative(y, x), 
   dy0.subs(-sols[0], -y))
Out[34]:
$\displaystyle \frac{d}{d x} y = – \frac{x}{y}$

2つ目の解についても同様に。

In [35]:
Eq(y, sols[1])
Out[35]:
$\displaystyle y = \sqrt{1 – x^{2}}$
In [36]:
dy1 = diff(sols[1], x)
Eq(Derivative(y, x), dy1)
Out[36]:
$\displaystyle \frac{d}{d x} y = – \frac{x}{\sqrt{1 – x^{2}}}$
In [37]:
Eq(Derivative(y, x), 
   dy1.subs(sols[1], y))
Out[37]:
$\displaystyle \frac{d}{d x} y = – \frac{x}{y}$

以上のことから，1つ目の解 sols[0] についても2つ目の解 sols[1] についても，

$$\frac{dy}{dx} = – \frac{x}{y}$$

となることがわかる。これは陰関数定理で求めた解と一致する。

例題

$f(x,y) = x^2 + y^2 -1 = 0$ を満たす陰関数 $y = y(x)$ の極大値と極小値を求めよ。

極値を求めるには $\displaystyle \frac{dy}{dx}$ と極大・極小（上に凸・下に凸）を判断するために $\displaystyle \frac{d^2y}{dx^2}$ が必要。

陰関数の1階微分

すでに，$\displaystyle \frac{dy}{dx} = – \frac{x}{y}$ は求めていて，変数 sol[0] に格納している。これを dy として使っていこう。

In [38]:
y = Function('y')
dy = sol[0]
Eq(Derivative(y(x), x), dy)
Out[38]:
$\displaystyle \frac{d}{d x} y{\left(x \right)} = – \frac{x}{y{\left(x \right)}}$
陰関数の2階微分

dy すなわち sol[0] をもう1階 $x$ で微分して $\displaystyle \frac{d^2y}{dx^2}$ を求める。

In [39]:
ddy = diff(dy, x)
Eq(Derivative(y(x), x, 2), ddy)
Out[39]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = \frac{x \frac{d}{d x} y{\left(x \right)}}{y^{2}{\left(x \right)}} – \frac{1}{y{\left(x \right)}}$

上式の右辺に dy つまり $\displaystyle \frac{dy(x)}{dx} = – \frac{x}{y(x)}$ を代入してやると…

In [40]:
ddy = ddy.subs(Derivative(y(x), x), dy)
Eq(Derivative(y(x), x, 2), ddy)
Out[40]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = – \frac{x^{2}}{y^{3}{\left(x \right)}} – \frac{1}{y{\left(x \right)}}$

極値となりそうな $x$ の値は

$\displaystyle \frac{dy}{dx} = – \frac{x}{y} = 0$ より $x = 0$.

そのときの $y$ の値は連立方程式 $x = 0, \ f(x, y)=0$ を解いて…

In [41]:
ans2 = solve([x, x**2 + y(x)**2 - 1], [x, y(x)])
ans2
Out[41]:
[(0, -1), (0, 1)]

ans2[0] つまり $(x, y) = (0, -1)$ のとき，$\displaystyle \frac{d^2y}{dx^2}$ つまり ddy の右辺の値を評価すると…

In [42]:
Eq(Derivative(y(x), x, 2),
   ddy.subs(x, ans2[0][0]).subs(y(0), ans2[0][1]))
Out[42]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = 1$

… となり，$\displaystyle \frac{d^2y}{dx^2} = 1 > 0$ だからここでは下に凸，つまり極小値（最小値）。

同様に ans2[1] つまり $(x, y) = (0, 1)$ の場合は…

In [43]:
Eq(Derivative(y(x), x, 2),
   ddy.subs(x, ans2[1][0]).subs(y(0), ans2[1][1]))
Out[43]:
$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)} = -1$

… となり，$\displaystyle \frac{d^2y}{dx^2} = -1 < 0$ だからここでは上に凸，つまり極大値（最大値）。

参考：陰関数のグラフ

SymPy Plotting Backends では $f(x, y) = 0$ という陰関数表示のままでグラフを描くことができる。

In [44]:
var('x y f')
f = x**2 + y**2 - 1

p1 = plot_implicit(f, "$f(x,y)=0$", 
                   (x, -1.2, 1.2), (y, -1.2, 1.5), 
                   legend = True,
                   size=(6,6),
                   aspect = "equal");

$y(x)$ の最大値，最小値の点もあわせてグラフに。（以前のバージョンでは plot_implicit() 以外の凡例が表示されないという不具合があったが，v2.1.0 で修正された。）

In [45]:
p2 = plot_list([0],[1], "y(x) の最大値",
               is_point = True, legend = True,
               xlim=(-1.2, 1.2), ylim=(-1.2, 1.5), 
               size=(6,6),
               aspect = "equal", show = False)

p3 = plot_list([0],[-1], "y(x) の最小値",
               is_point = True, legend = True,
               xlim=(-1.2, 1.2), ylim=(-1.2, 1.5), 
               size=(6,6),
               aspect = "equal", show = False)

p = p1+p2+p3
p.show();