\usepackagecancel

Return to SymPy 演習

3. SymPy で微分・積分

1変数関数,特に(双曲線関数・逆双曲線関数を含む)初等関数の微分・積分。また,2変数関数に関する偏微分と2重責分,および2重積分に関連したガウス積分について。

SymPy の import

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

1変数関数の微分 diff()

1階微分

SymPy で1変数関数 y=f(x) の微分には以下のように diff() を使います。

dfdx=ddxf(x)= diff(f(x), x)= Derivative(f(x), x).doit()

高階微分

2階以上の高階微分は

dnfdxn=dndxnf(x)= diff(f(x), x, n)

参考:Derivative()diff()

Derivative() は微分の表示をするだけで実際に微分は行いません。微分を実行するには,Derivative().doit() のように .doit() を追加します。

diff() を使えば即微分を実行します。

例:

y=x3+5x+2x で微分します。とっとと答えだけを表示させるには…

In [2]:
diff(x**3 + 5*x + 2, x)
Out[2]:
3x2+5

何を計算させるのかを表示させてから答えを表示させるには…

In [3]:
y = x**3 + 5*x + 2

Eq(Derivative(y, x), 
   Derivative(y, x).doit())
Out[3]:
ddx(x3+5x+2)=3x2+5

y=x3+5x+2x で2階微分します。

In [4]:
Eq(Derivative(y, x, 2), diff(y, x, 2))
Out[4]:
d2dx2(x3+5x+2)=6x

べき関数,平方根の微分

べき関数 y=xp の微分

In [5]:
diff(x**p, x)
Out[5]:
pxpx

直前の結果を _ で参照して簡単化 simplify() してみます。

In [6]:
simplify(_)
Out[6]:
pxp1

ということで,

In [7]:
Eq(Derivative(x**p, x), diff(x**p, x).simplify())
Out[7]:
xxp=pxp1

偏微分記号が出てくるのは,x の他に p も「変数」とみなされているから。

特に平方根 x の微分は…

ddxx= diff(sqrt(x), x)

In [8]:
Eq(diff(sqrt(x), x, evaluate = False), 
   diff(sqrt(x), x))
Out[8]:
ddxx=12x

指数関数の微分

ネイピア数 e を底とした指数関数 ex= exp(x) の微分

ddxex= diff(exp(x), x)

In [9]:
diff(exp(x), x)
Out[9]:
ex

対数関数の微分

自然対数 log(x) の微分

ddxlogx= diff(log(x), x)

In [10]:
diff(log(x), x)
Out[10]:
1x

三角関数の微分

sinx, cosx, tanx の微分

In [11]:
Eq(Derivative(sin(x), x), diff(sin(x), x))
Out[11]:
ddxsin(x)=cos(x)
In [12]:
Eq(Derivative(cos(x), x), diff(cos(x), x))
Out[12]:
ddxcos(x)=sin(x)
In [13]:
Eq(Derivative(tan(x), x), diff(tan(x), x))
Out[13]:
ddxtan(x)=tan2(x)+1

右辺を簡単化します。
diff(tan(x), x).simplify()simplify(diff(tan(x), x)) と同じ意味です。

In [14]:
Eq(Derivative(tan(x), x), diff(tan(x), x).simplify())
Out[14]:
ddxtan(x)=1cos2(x)

逆三角関数の微分

ddxsin1x=ddxarcsinx= diff(asin(x), x)

In [15]:
Eq(Derivative(asin(x), x), diff(asin(x), x))
Out[15]:
ddxasin(x)=11x2

○練習:arccosx, arctanx の微分

ddxcos1x=ddxarccosx= diff(acos(x), x)

ddxtan1x=ddxarctanx= diff(atan(x), x) についても同様にやってみてください。

In [ ]:
 
In [ ]:
 

双曲線関数の微分

ddxcoshx= diff(cosh(x), x)

In [16]:
Eq(Derivative(cosh(x), x), diff(cosh(x), x))
Out[16]:
ddxcosh(x)=sinh(x)

○練習:sinhx, tanhx の微分

ddxsinhx= diff(sinh(x), x)

ddxtanhx= diff(tanh(x), x) についても同様に。

In [ ]:
 
In [ ]:
 

逆双曲線関数の微分

ddxcosh1x=ddxarcosh x= diff(acosh(x), x) =1x21

In [17]:
Eq(Derivative(acosh(x), x), diff(acosh(x), x))
Out[17]:
ddxacosh(x)=1x1x+1

本稿執筆時点では,SymPy は 1x21 ではなく,あえて 1x1x+1 と返すようです。(私見ですが,これが後々,いくつかの残念な結果を引き起こします。)

○練習:sinh1x, tanh1x の微分

ddxsinh1xddxtanh1x についても同様に微分してみて。

In [ ]:
 
In [ ]:
 

1変数関数の積分 integrate()

不定積分

SymPy で1変数関数 y=f(x) の積分には integrate() を使います。

f(x)dx= integrate(f(x), x)= Integral(f(x), x).doit()

不定積分の際,積分定数が省略された答えが出力されます。

定積分

定積分の際は以下のように積分範囲を (x, a, b) のように指定します。

abf(x)dx= integrate(f(x), (x, a, b))

参考:Integral()integrate()

Integral() は積分の表示をするだけで実際に積分は行いません。積分を実行するには,Integral().doit() のように .doit() を追加します。

integrate() を使えば即積分を実行します。

べき関数の不定積分

xpdx= integrate(x**p, x)

In [18]:
Eq(Integral(x**p, x), integrate(x**p, x))
Out[18]:
xpdx={xp+1p+1forp1log(x)otherwise

p=1 のときは,以下の「対数関数に関連した不定積分」も参照:

指数関数の不定積分

exdx= integrate(exp(x), x)

In [19]:
Eq(Integral(exp(x), x), integrate(exp(x), x))
Out[19]:
exdx=ex

対数関数に関連した不定積分

1xdx=log|x|+C となるはずですが…

In [20]:
Eq(Integral(1/x, x), integrate(1/x, x))
Out[20]:
1xdx=log(x)

SymPy の積分 integrate() では log の中身の絶対値を省略する傾向にある。

○練習:logx の不定積分

logxdx= integrate(log(x), x) =

In [ ]:
 

上記の結果は,いかにも部分積分をして答えを出しているようですね。

三角関数の不定積分

In [21]:
Eq(Integral(cos(x), x), 
   Integral(cos(x), x).doit())
Out[21]:
cos(x)dx=sin(x)
In [22]:
Eq(Integral(sin(x), x), 
   Integral(sin(x), x).doit())
Out[22]:
sin(x)dx=cos(x)
In [23]:
Eq(Integral(1/cos(x)**2, x), 
   integrate(1/cos(x)**2, x).simplify())
Out[23]:
1cos2(x)dx=tan(x)

○練習:tanx の不定積分

では,tanxdx は?

In [ ]:
 

逆三角関数に関連した不定積分

ddxarcsinx=11x2 だったので,

11x2dx=arcsinx+C となるはず。

In [24]:
Eq(Integral(1/(sqrt(1 - x**2)), x), 
   integrate(1/(sqrt(1 - x**2)), x))
Out[24]:
11x2dx=asin(x)

ddxarccosx=11x2 だったので,

11x2dx=arccosx+C となってもいいはずだが…

In [25]:
Eq(Integral(-1/(sqrt(1 - x**2)), x), 
   integrate(-1/(sqrt(1 - x**2)), x))
Out[25]:
(11x2)dx=asin(x)

積分の答えが arcsinx と出ても,慌てず騒がず,直近の出力 _.rewrite(acos) つまり arccosx を使って書き直してみると…

In [26]:
_.rewrite(acos)
Out[26]:
(11x2)dx=acos(x)π2

π2 は積分定数 C に組み込まれる定数なので,結果,

11x2dx=arccosx+C となってもいいことになる。

ddxarctanx=11+x2 だったので,

11+x2dx=arctanx+C となるはず。

In [27]:
integrate(1/(1 + x**2), x)
Out[27]:
atan(x)

○練習:逆三角関数の不定積分

arcsinxdx, arccosxdx, arctanxdx についても計算してみよ。

In [ ]:
 
In [ ]:
 
In [ ]:
 

双曲線関数の不定積分

In [28]:
Eq(Integral(cosh(x), x), integrate(cosh(x), x))
Out[28]:
cosh(x)dx=sinh(x)
In [29]:
Eq(Integral(sinh(x), x), integrate(sinh(x), x))
Out[29]:
sinh(x)dx=cosh(x)
In [30]:
Eq(Integral(1/cosh(x)**2, x), integrate(1/cosh(x)**2, x).simplify())
Out[30]:
1cosh2(x)dx=tanh(x)

○練習:tanhx の不定積分

では tanhxdx は?

In [ ]:
 

逆双曲線関数に関連した不定積分

1. 1x21dx

ddxcosh1x=ddxarcosh x= diff(acosh(x), x) =1x21 なので

1x21dx=cosh1x となるはずだが…

In [31]:
Eq(Integral(1/sqrt(x**2 - 1), x), 
   integrate(1/sqrt(x**2 - 1), x))
Out[31]:
1x21dx=log(x+x21)

cosh1xlog(x+x21) と書き直してはくれるが…

In [32]:
acosh(x).rewrite(log).factor()
Out[32]:
log(x+x1x+1)

log(x+x21)cosh1x とは書き直してくれない。

In [33]:
log(x + sqrt(x**2 - 1)).rewrite(log)
Out[33]:
log(x+x21)

2. 11+x2dx

ddxsinh1x=11+x2 なので…

In [34]:
Eq(Integral(1/sqrt(1 + x**2), x), 
   integrate(1/sqrt(1 + x**2), x))
Out[34]:
1x2+1dx=asinh(x)

3. 11x2dx

ddxtanh1x=11x2 なので

11x2dx=tanh1x+C となるはずだが…

In [35]:
Eq(Integral(1/(1-x**2), x), integrate(1/(1-x**2), x))
Out[35]:
11x2dx=log(x1)2+log(x+1)2

4. sinh1xdx

In [36]:
Eq(Integral(asinh(x), x), integrate(asinh(x), x))
Out[36]:
asinh(x)dx=xasinh(x)x2+1

5. tanh1xdx

In [37]:
Eq(Integral(atanh(x), x), integrate(atanh(x), x))
Out[37]:
atanh(x)dx=xatanh(x)+log(x+1)atanh(x)

6. cosh1xdx

In [38]:
integrate(acosh(x), x)
Out[38]:
acosh(x)dx

上記のように,本稿執筆時点では,SymPy は cosh1x の積分ができない。これは部分積分で簡単に求まるから,以下のように人力で計算してみよう。

cosh1xdx=(x)cosh1xdx=xcosh1xx(cosh1x)dx=

SymPy がこの積分ができないのは,どうも以下のように…

In [39]:
Eq(Derivative(acosh(x), x), diff(acosh(x), x))
Out[39]:
ddxacosh(x)=1x1x+1

ddxcosh1x1x21 ではなく 1x1x+1 とすることにあるのでは?と推測する。

ためしに,xx1x+1dx を計算してみると…

In [40]:
integrate(x/(sqrt(x-1)*sqrt(x+1)), x)
Out[40]:
G6,66,2(14,140,0,12,112,14,0,14,12,0|1x2)4π32+iG6,62,6(1,34,12,14,0,134,141,12,12,0|e2iπx2)4π32

なんか意味不明の答えが出てくるであろう。

一方,xx21dx を計算してみると…

In [41]:
Eq(Integral(x/sqrt(x**2-1), x), 
   integrate(x/sqrt(x**2-1), x))
Out[41]:
xx21dx=x21

以上のことから

cosh1xdx=xcosh1xx(cosh1x)dx=xcosh1xxx21dx=

○練習:cosh1x の不定積分

部分積分を使って,あらためて cosh1xdx を求めよ。

In [ ]:

偏微分:多変数(主に2変数)関数の微分

主に2変数関数 f(x,y) の偏微分についてまとめておきます。

1階偏微分 diff()

SymPy での偏微分は(1変数関数の微分である「常微分」と同様の書き方で)以下のように書きます。

xf(x,y)= diff(f(x, y), x)

yf(x,y)= diff(f(x, y), y)

In [42]:
# 変数 x, y を初期化
var('x y')

# f を未定義関数として宣言
f = Function('f')
In [43]:
diff(f(x, y), x)
Out[43]:
xf(x,y)
In [44]:
diff(f(x, y), y)
Out[44]:
yf(x,y)

○練習:偏導関数を求める

関数 z=x34x2y+xy+3y2 の偏導関数を求めよ。

「偏導関数を求めよ」とは,x,y の2変数関数である z について zx
zy を計算して求めよ,ということ。

In [ ]:
In [ ]:
 

高階偏微分 diff(f(x, y), x, m, y, n)

2階以上の高階偏微分は以下のように書きます。

2x2f(x,y)= diff(f(x, y), x, 2)

y(xf(x,y))=2yxf(x,y)= diff(f(x, y), x, y)

nyn(mxmf(x,y))= diff(f(x, y), x, m, y, n)

○練習:高階導関数の計算

1.  2変数関数 z=tan1yx=arctanyx について,以下の高階偏導関数を計算し,簡単化した表示で表せ。

(1). 2zyx,(2). 2zxy,(3). 2zx2,(4). 2zy2

In [ ]:
 
In [ ]:
 
In [ ]:
 

2. z=logx2+y2 のとき,次の値を求めよ。

2zx2+2zy2=?

In [ ]:

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

2fxy=2fyx

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

以下では,比較演算子 == を使って左辺と右辺が等しいかを確認しています。等しければ True,等しくなければ False を返します。

In [45]:
diff(f(x, y), x, y) == diff(f(x, y), y, x)
Out[45]:
True

以下では,左辺から右辺を引いてゼロになるかどうかで,確認しています。

In [46]:
Eq(Derivative(f(x, y), x, y)-Derivative(f(x, y), y, x),
   diff(f(x, y), x, y)-diff(f(x, y), y, x))
Out[46]:
2yxf(x,y)2xyf(x,y)=0

2重積分:2変数関数の積分

多変数関数の積分である多重積分のうち,2変数関数の積分である2重積分についてまとめておきます。

累次積分

2重積分の計算は,基本的に累次積分でやります。

Df(x,y)dxdy の領域 D が何と何で囲まれているかを明らかにして積分する。

ケース 1

領域 Dx1xx2, y1yy2 である場合:

Df(x,y)dxdy=x1x2{y1y2f(x,y)dy}dx

または

Df(x,y)dxdy=y1y2{x1x2f(x,y)dx}dy

ケース 2

領域 Dx1xx2, y1(x)yy2(x) である場合:

Df(x,y)dxdy=x1x2{y1(x)y2(x)f(x,y)dy}dx

ケース 3

領域 Dx1(y)xx2(y), y1yy2 である場合:

Df(x,y)dxdy=y1y2{x1(y)x2(y)f(x,y)dx}dy

SymPy では2重積分は以下のように:

y1y2{x1x2f(x,y)dx}dy=
integrate(f(x, y), (x, x1, x2), (y, y1, y2))

○練習:累次積分による2重積分

I=Ddxdy,D:x2+y21

累次積分の形にする。

領域 D の定義式を y について解き,

1x2y1x2

平方根の中がゼロ以上であることから

1<x<1

被積分関数は f(x,y)=1

ケース 2 の累次積分の形:
f(x,y)=1x1=1x2=1y1(x)=1x2y2(x)=1x2

Df(x,y)dxdy=x1x2{y1(x)y2(x)f(x,y)dy}dx=

In [47]:
f = 1
x1 = -1
x2 = 1
y1 = -sqrt(1-x**2)
y2 = sqrt(1-x**2)

Integral(f, (y, y1, y2), (x, x1, x2))
Out[47]:
111x21x21dydx

確認のために,段階ごとに計算してみます。

まず,1x21x2fdy= integrate(f, (y, y1, y2))

In [48]:
Eq(Integral(f, (y, y1, y2)), 
   integrate(f, (y, y1, y2)))
Out[48]:
1x21x21dy=21x2

次に,1121x2dx を計算します。

人力でこの積分を計算するには,以下の不定積分の項を参照のこと。

人力だと一手間も二手間もかかりますが,SymPy だと簡単に積分してくれます。

積分を続けて最終的な答えを求めなさい。

In [ ]:
In [ ]:

ガウス積分の準備としての2重積分

S=ex2y2dxdy

In [49]:
Eq(Integral(exp(-x**2-y**2), (x, -oo, oo), (y, -oo, oo)),
   integrate(exp(-x**2-y**2), (x, -oo, oo), (y, -oo, oo)))
Out[49]:
ex2y2dxdy=π

SymPy では上記のように苦もなく計算しているが,人力で積分するためには以下のように極座標に変換して積分することになる。

極座標に変換して積分

極座標に変換して計算する。

デカルト座標 x,y と2次元極座標 r,θ との間の関係は
x=rcosθy=rsinθ
であり,微小面積要素 dxdy は以下のように変換される。

dxdy=(x,y)(r,θ)drdθ

ここで,ここで,(x,y)(r,θ)|xrxθyryθ|=xryθyrxθ をヤコビアンという。

SymPy ではヤコビ行列を計算する機能があるが,簡単なのでここでは以下のように計算してみる。

In [50]:
x = r*cos(theta)
y = r*sin(theta)

jaco = diff(x, r)*diff(y,theta) - diff(y,r)*diff(x,theta)
simplify(jaco)
Out[50]:
r

ヤコビアンが r となることがわかった。したがって極座標への変換によって

x2+y2=r2ex2y2=e(x2+y2)=er2dxdy=rdrdθD:{<x<,<y<}=D:{0r<,0θ2π}

となり,

  S=ex2y2dxdy=02π{0er2rdr}dθ=02πdθ0er2rdr=2π0er2rdr

であるから,この問題は 0er2rdr ができればよい。これも SymPy は何なく積分してしまう。

In [51]:
Eq(Integral(r*exp(-r**2), (r, 0, oo)), 
   integrate(r*exp(-r**2), (r, 0, oo)))
Out[51]:
0rer2dr=12

人力で積分する際はきっと r2=t とおいて置換積分するんですよねぇ。

参考:敢えて置換積分 .transform()

せっかくだから,これを敢えて置換積分でやってみる。

まず,この積分を int1 として定義。

In [52]:
int1 = Integral(r*exp(-r**2), (r, 0, oo))
int1
Out[52]:
0rer2dr

r2t という変数に置換する。

In [53]:
int2 = int1.transform(r**2, t)
int2
Out[53]:
0et2dt

この形なら暗算でも積分できるが,念のため .doit() で実際に積分を実行させる。

In [54]:
Eq(int2, int2.doit())
Out[54]:
0et2dt=12

したがって

S=2π0er2rdr=2π×12=π

ガウス積分

I=ex2dx

とすると,上で

S=ex2y2dxdy=I2=π

ということがわかったので,

I=ex2dx=S=π

ということになる。SymPy がガウス積分 I を知っていることを確認する。

In [55]:
# 変数 x の初期化
var('x')

Eq(Integral(exp(-x**2), (x, -oo, oo)),
   Integral(exp(-x**2), (x, -oo, oo)).doit())
Out[55]:
ex2dx=π

被積分関数 ex2 は偶関数であるから,積分範囲を 0 からにすると答えは半分になるはずで,これも確認。

In [56]:
Eq(Integral(exp(-x**2), (x, 0, oo)),
   Integral(exp(-x**2), (x, 0, oo)).doit())
Out[56]:
0ex2dx=π2

ガウス積分は,積分の上限が のときにだけ解析的に積分できる例。不定積分だと積分できない(結果が初等関数で表せない)。
… のだが,実際に SymPy でやってみると…

In [57]:
Eq(Integral(exp(-x**2), x),
   Integral(exp(-x**2), x).doit())
Out[57]:
ex2dx=πerf(x)2

参考:誤差関数

erf(x) は「誤差関数」と呼ばれる特殊関数(つまり初等関数では表せないということ)で,これが定義:

erf(x)2π0xet2dt

○練習:誤差関数の極限値

誤差関数 erf(x)x でガウス積分で書けるようになるので,

limxerf(x)=2π0et2dt=2ππ2=1

となるはずである。これを SymPy で確認してください。

In [ ]: