関数，微積分，級数展開

Sympyの基礎

Sympyは数式や変数を記号のまま代数的に計算できるPythonのライブラリである．数理解析で学ぶ微分方程式・フーリエ解析をSympyで理解するための最小限の使い方をここでまとめる．

まず，変数（variable） とは，任意の範囲で任意の値を取ることができる記号（文字）である．具体的には，\(y=2x+1\)という一次関数があったとき，変数は \(x\) と \(y\) である．

Pythonにおける変数と計算の扱いを確認する．

x = 1
y = 2 * x + 1

print(f'x={x}, y={y}')

x=1, y=3

\(y=2x+1\) という式に \(x=1\) を代入したときの結果を計算できた．

Caution

プログラミングにおける変数と数式の変数は意味が異なる点に注意されたい．以下のようにプログラミングでは，変数に値を代入し，値をメモリ上に確保している．データの型に応じて任意の値を代入可能であるが，\(y=2x+1\)という数式に現れる変数とは違い，すでに値を代入して式として扱っている点が異なる．

例えば，\([-1,1]\) の範囲で，\(y=2x+1\) を描画したいとき，この範囲の値を各代入する必要がある．

numpy を使って次のようにプログラミングする．

import numpy as np
x = np.linspace(-1, 1, 100)
y = 2 * x + 1

Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.

上記のセルでは，np.linspace を使って，\([-1,1]\) の範囲で \(100\) 個の数列を作成し，この数列を y の式に代入している．

そのため，以下のように .shape で確認すると，x と y は複数の値から構成されていることがわかる．

print(x.shape)
print(y.shape)

(100,)
(100,)

matplotlib を使って可視化しよう．

import matplotlib.pyplot as plt

plt.plot(x, y, label='y = 2x + 1')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(True)
plt.legend()
plt.show()

本講義で扱うSympyでは，代入をすることなく，文字として変数を扱ったプログラミングが可能である．つまり，これまでの数学の授業で扱ってきたように数式と変数を定義できる．まずはSympyから変数を扱うためのsymbolsをimportする．

from sympy import symbols

Sympyでは以下のように変数を定義する．printすると文字式の変数として \(x\) が定義されていることがわかる．変数 \(x\) には上記の\(x=1\)のような具体的な値が代入されていないことを確認されたい．

x = symbols('x')

x

\[\displaystyle x\]

\(y=2x+1\) を作成する．

y = 2 * x + 1
y

\[\displaystyle 2 x + 1\]

式を作ることができた．これを可視化する．sympyでは便利なplot関数が用意されており，定義した関数を容易にプロットすることができる．詳細はsympyのplottingモジュールを参考にされたい．

from sympy.plotting import plot
plot(y, (x, -1, 1))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x1248928b0>

これまでのように変数 \(x\) に具体的な値を代入したい場合は，メソッド 式.subs(変数, 代入したい値)を利用する．

x.subs(x, 1)

\[\displaystyle 1\]

y は x の式なので次のように代入もできる．

y.subs(x, 1)

\[\displaystyle 3\]

\(2 \times 1 + 1 = 3\) と正しく計算されていることがわかる．

関数

関数（function） とは，ある変数 \(x\) について，一つの変数 \(y\) に対応させるものである．

\[ y = f(x) \]

ここで \(f(\cdot)\) を関数と呼び，変数 \(x\) が取りうる範囲を 定義域（domain） ，変数 \(y\) が取りうる値を 値域（range） と言う．

具体的には以下のようなものを関数と呼ぶ．

\[\begin{split} \begin{align} &y = 2x + 1, \\ &y = x^2 - 1, \\ &x^2 + y^2 = 1 \end{align} \end{split}\]

他にもこれまで様々な関数と出会ってきたはずである．また前述したような関数は変数を一つだけ取りうるため一変数関数と呼ばれ，多変数を受け取る多変数関数も考えることができる（後述）．

以降で，Sympyを使って変数を定義して，関数を定義することを試みる．

一変数関数

では，定義した\(x\) を使って一変数関数を定義する．一変数関数とは，関数 \(f(\cdot)\) において取りうる変数の数が一つだけの関数を示す．

\[ y = f(x) \]

これは先ほど実装したが，もう一度定義しておく．

y = 2 * x + 1

y

\[\displaystyle 2 x + 1\]

等式 \(=\) で定義したい場合，Eq クラスをインポートして次のように定義する．

from sympy import Eq

y = symbols('y')
eq = Eq(y, 2 * x + 1)
eq

\[\displaystyle y = 2 x + 1\]

y を symbol として定義する必要がある点に注意してほしい．微分方程式を扱う際に，Eq を利用するので覚えておこう．

多変数関数

2以上の正の実数 \(n\) に対して，\(n\) 個の変数 \(x_1, ..., x_n\) を持つ関数　\(f(x_1, x_2, ..., x_n)\) を多変数関数と呼ぶ．

\[ y = f(x_1, x_2, ..., x_n) \]

例えば，以下のような2変数 \(x, y\) を持つ多変数関数 \(f(x, y)\) を考えることができる．

1. \( f(x_1, x_2) = x_{1}^2 + x_{2}^2 \)

2. \( f(x_1, x_2) = \cos(x_1) + \sin(x_2) \)

3. \( f(x_1, x_2, x_3) = x_{1}^3 + x_{2}^2 + x_{3} + 1 \)

では，数式 \(y=x_{1}^2 + x_{2}^2\) を実装する．まずは変数 \(x_1, x_2\) を定義する．

x_1 = symbols('x_1')
x_2 = symbols('x_2')

なお，symbol は次のように複数同時に定義できる

x_1, x_2 = symbols('x_1 x_2')

続いて，二変数の入力を持つ多変数関数 \(f(x_1, x_2)\) を定義する．

y = x_1**2 + x_2**2

y

\[\displaystyle x_{1}^{2} + x_{2}^{2}\]

定義することができた．一変数関数と同様に代入を行うが，多変数なのでsubsメソッドには辞書形式の引数として与える必要がある．

y.subs({x_1: 1, x_2: 2})

\[\displaystyle 5\]

次のようにも代入できるが，可読性の面から上記の代入方法を推奨する．

y.subs(x_1, 1).subs(x_2, 2)

\[\displaystyle 5\]

関数のプロット

関数を扱う際，関数がどのような形を持つかを可視化することは理解の助けとなる．一変数関数と多変数関数を可視化しよう．

一変数関数のプロット

関数 \(y=x^2+2x+1\) を定義して，定義域を \([-3, 3]\) でプロットする．

from sympy.plotting import plot
x = symbols('x')
y = x**2 + 2*x + 1
plot(y, (x, -3, 3))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x128501d30>

二変数関数のプロット

関数 \(y=x_1^2+x_2^2+1\) を定義して，\(x_1\) と \(x_2\) のそれぞれの定義域を \([-3, 3]\) でプロットする．

from sympy.plotting import plot3d
x_1 = symbols('x_1')
x_2 = symbols('x_2')
y = x_1**2 + x_2**2 + 1
plot3d(y, (x_1, -3, 3), (x_2, -3, 3))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x1249e23a0>

Function クラス

微分方程式では，変数 t に依存する未知関数 x(t) そのものを定義したい場面がある．その場合，Function クラスを使う．

from sympy import Function

t = symbols('t')
x = Function('x')
x_t = x(t)
x_t

\[\displaystyle x{\left(t \right)}\]

次のようにも実装できる．

x_t = Function('x')(t)
x_t

\[\displaystyle x{\left(t \right)}\]

微積分

本講義で扱う微分方程式に 微分 という言葉が含まれていることからもわかるように，微分は関数の性質を知る上で非常に重要な役割を持つ．まずは一変数関数から微分の意味を復習する．

微分

これまでの中学〜大学の講義で微分とは関数の傾き・変化量であると習ったはずである．この関数の傾きを知るために，点 \(x_1\) で関数 \(f\) に接する直線 \(g\) を考える．この直線 \(g\) は接線であった．接線の傾きを知るために，ある二点 \(x_1,x_2\)を通る直線の傾きを考える．この傾き \(a\) は以下の式に示すように増加量で計算することができる．

\[ a=\frac{f\left(x_{2}\right)-f\left(x_{1}\right)}{x_{2}-x_{1}} \]

\(x_1\)での傾きを知るために，\(x_2\)を\(x_1\)に限りなく近づける．\(x_1\)を\(x\)，\(x_2-x_1\)を\(h\)として極限を使って表すと，

\[ \lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h} \]

となる．この極限の値は，\(f\)の 導関数（derivative） と呼び，\(f'(x)\) または \(\frac{df}{dx}\)と書く．この導関数を求めることを 微分 するという．

また\(f^{(n)}\)は\(n\)階の導関数であり，微分方程式に含まれる導関数の最大の次数 \(n\) を微分方程式の階数と呼ぶ．

sympyではdiffメソッドを利用すると容易に微分が計算できる．関数 \(y=x^2+2x+1\) の一次微分 \(f'(x)\) を計算する．

from sympy import diff
x = symbols('x')
f = x**2 + 2*x + 1
f_x = diff(f, x)

f_x

\[\displaystyle 2 x + 2\]

2階微分は次のように計算できる

f_xx = diff(f, x, 2)
f_xx

\[\displaystyle 2\]

微分した結果が正しく計算されたことを手計算と併せて確認されたい．

点 \(x'\) で関数 \(f\) に接する直線（接線） \(g\) を求める．微分して接線の傾きが\(a\)のとき，\(g\)は以下の式で計算できることを思い出して欲しい．

\[ g = f(x') + a (x - x') \]

x_ = 1
g = f.subs({x:x_}) + f_x.subs({x:x_}) * (x - x_)

g

\[\displaystyle 4 x\]

計算された接線を可視化してみる．

fig = plot(f, g, (x, -10, 10), xlim=(-5, 5), ylim=(-1, 10), show=False)
fig[0].line_color = 'b'
fig[1].line_color = 'r'
fig.show()

色は fig[0].line_color = 'b' のように指定できる．

偏微分

多変数関数において，ある入力 \(x_i\) で微分することを 偏微分 (partial differentiation) と呼ぶ．変数 \(x,y\) を入力として受け取る多変数関数 \(f(x,y)\) の偏微分は以下の式で与えられ，偏微分して得られる関数を 偏導関数 (partial derivative) と呼ぶ．

\[\begin{split} \begin{align} \frac{\partial f}{\partial x}(x, y)&=\lim_{h \rightarrow 0} \frac{f(x+h, y)-f(x, y)}{h} \\ \frac{\partial f}{\partial y}(x, y)&=\lim_{h \rightarrow 0} \frac{f(x, y+h)-f(x, y)}{h} \end{align} \end{split}\]

偏微分を使うと，ある点で局面に接する平面，接平面(tangent plane) を書くことができる．そのために，勾配 (gradient) \(\nabla_f\) を定義する．

\[ \nabla_f = (\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}) \]

偏微分も同様にdiffメソッドで第二引数に微分対象の変数を指定すれば実行できる．

x = symbols('x')
y = symbols('y')
f = x**2 + y**2 + 1

f_x = diff(f, x)
f_y = diff(f, y)

f_x

\[\displaystyle 2 x\]

f_y

\[\displaystyle 2 y\]

微分した結果が正しく計算されたことを手計算と併せて確認されたい．

接平面を求める．関数 \(f\) の点(a, b)における接平面の方程式は以下で与えられる．

\[ g(x, y) = f_x(a, b) (x - a) + f_y(a, b) (y - b) + f(a, b) \]

p = (1, 1)
a, b = p
g = f_x.subs({x:a, y:b}) * (x - a) + f_y.subs({x:a, y:b}) * (y - b) +  f.subs({x:a, y:b})
g

\[\displaystyle 2 x + 2 y - 1\]

plot3d(f, g, (x, -5, 5), (y, -5, 5))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x1295587f0>

Derivative クラス

微分は diff() メソッドを呼び出すと計算することができる．一方，微分方程式をプログラミングする際，導関数・偏導関数そのものを定義したい場合がある．このようなときに利用するのが Derivative クラスである．次のように利用する．

from sympy import Derivative

f = x**2 + 2*x + 1
df = Derivative(f, x)
df

\[\displaystyle \frac{d}{d x} \left(x^{2} + 2 x + 1\right)\]

微分を評価していない導関数を定義できた．ここから微分したい場合は doit() を使う．

df.doit()

\[\displaystyle 2 x + 2\]

前述の結果が確認できたかと思う．Eq を使うと等式で結べる．

Eq(df, diff(f, x))

\[\displaystyle \frac{d}{d x} \left(x^{2} + 2 x + 1\right) = 2 x + 2\]

Function と Eq クラスと組み合わせると微分方程式を定義できる．

from sympy import Eq, Function, Derivative, symbols

t = symbols('t')
x = Function('x')(t)
d_x = Derivative(x, t)
ode = Eq(d_x, -x)
ode

\[\displaystyle \frac{d}{d t} x{\left(t \right)} = - x{\left(t \right)}\]

積分

積分を利用することで関数の面積や曲線の長さを計算することができた．積分には積分範囲や被積分関数 に応じていくつかの種類があるが，ここでは数理解析の講義で必要となる不定積分，定積分，広義積分を復習し，実装する．

以下の一変数の二次関数を扱う．関数の定義と変数 \(x\) に関する微分を実行する

x = symbols('x')
f = x**2 + 2*x + 1
f

\[\displaystyle x^{2} + 2 x + 1\]

d_x = diff(f, x)
d_x

\[\displaystyle 2 x + 2\]

定積分

定積分は与えられた被積分関数 \(f(x)\) について，積分区間 \([a,b]\) でその関数の面積を求める方法であった．以下のように，積分変数 \(dx\) と積分 \(\int\) を使って定積分を書くことができる．

\[ F(x) = \int_{a}^{b} f(x) dx \]

sympyではintegrate関数を使うことで計算できる．被積分関数として微分された式を第一引数に，積分変数と範囲を第二引数として与える．

from sympy import integrate
f_ = integrate(d_x, (x, 0, 5))
f_

\[\displaystyle 35\]

不定積分

不定積分とは微分すると消える定数 \(C\) を積分定数として，元の被積分関数のすべての原始関数を計算する積分である．そのため，複数の解を持つので注意されたい．

\[ F(x) = \int f(x) dx + C \]

また不定積分は本講義で扱う微分方程式の解を求める際にも必要となるため復習されたい．sympyでは同様にintegrate関数を使って計算できる．

f_ = integrate(d_x, x)
f_

\[\displaystyle x^{2} + 2 x\]

広義積分

重要な積分として，広義積分がある．広義積分は積分区間が無限大や発散する点を含む場合に用いられる積分であり，フーリエ変換の導出の際に登場する以下のような積分である．

\[\begin{split} \begin{align} F(x) &= \int_{-\infty}^{\infty} f(x) dx, \\ F(x) &= \int_{a}^{\infty} f(x) dx, \\ F(x) &= \int_{-\infty}^{b} f(x) dx \end{align} \end{split}\]

ooは無限大を表す．

from sympy import oo
f_ = integrate(d_x, (x, 0, oo))
f_

\[\displaystyle \infty\]

from sympy import exp
f_ = integrate(exp(-x), (x, 0, oo))
f_

\[\displaystyle 1\]

Integrate クラス

Derivative と同様に，積分そのものを定義する際には Integral を利用する．

from sympy import Integral

f_ = Integral(d_x, x)
f_

\[\displaystyle \int \left(2 x + 2\right)\, dx\]

.doit() を使うと，積分を評価できる．

f_.doit()

\[\displaystyle x^{2} + 2 x\]

この機能はフーリエ変換で利用する．

補足：関数の連続性

関数が連続であるとは，極限 \(\lim_{x \rightarrow a} f(x)\) が　\(f(a)\) に一致すること，つまり，

\[ \lim _{x \rightarrow a} f(x) = f(a) \]

が成立するときに，関数 \(f(x)\) は点　\(a\) で連続であると言う．

\(a\) を1に近づけたとき，\(\lim_{x \rightarrow a} f(x) = f(a)\) が成立するか検証する

from sympy import limit
x = symbols('x')
f = x**2

a = 1.
g = limit(f, x, a)
g

\[\displaystyle 1.0\]

f.subs({x: a})

\[\displaystyle 1.0\]

級数展開

級数展開はある関数 \(f\) を関数の組み合わせで近似する手法である．本講義ではフーリエ解析の基礎として三角関数によるフーリエ級数展開を紹介する．その理解のために，関数を無限級数で近似するテイラー展開について復習する．

テイラー展開

テイラー展開は，関数を無限級数で近似する手法であり以下のように関数 \(f(x)\) をある点 \(a\) 周りで近似する．ただし，関数が無限回微分可能である必要性がある．

\[ f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)(x-a)^2}{2!} + \frac{f'''(a)(x-a)^3}{3!} + \cdots \]

ここで，\(f'\)，\(f''\)，\(f'''\)などはそれぞれ関数\(f\)の1次，2次，3次の導関数を表す．

テイラー展開はsympyでseries関数として実装されており，第一引数には関数を，第二引数に変数を，第三引数にどの点の周りで近似するかを，第四引数に展開の項数を与える．以下は\(f(x) = \sin(x)\)のテイラー展開の実装である（今回は \(x=0\) 周りでの展開なので正確にはマクローリン展開である）．

from sympy import sin, series
x = symbols('x')
f = sin(x)
f_taylor = series(f, x, 0, 7)

f_taylor

\[\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} + O\left(x^{7}\right)\]

ここで登場する removeO() はグラフプロットのために余剰項を削除するメソッドである．

plot(f, f_taylor.removeO(), (x, -3, 3))

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x129794c70>

補足：方程式

方程式を解く方法を紹介する．まず方程式は Eq を使って次のように定義できた．

eq = Eq(x**2 + 5*x + 6, 0)
eq

\[\displaystyle x^{2} + 5 x + 6 = 0\]

方程式の左辺と右辺は .lhs と .rhs で取得できる．

eq.lhs

\[\displaystyle x^{2} + 5 x + 6\]

eq.rhs

\[\displaystyle 0\]

この二次方程式は solve を使うと，方程式 eq を変数 x に関して解くことができる．

from sympy import solve
solutions = solve(eq, x)
solutions

[-3, -2]

Sympyは因数分解も factor でできるのでこちらからも確認する．

from sympy import factor
factored_eq = factor(eq.lhs)
factored_eq

\[\displaystyle \left(x + 2\right) \left(x + 3\right)\]

方程式が解けていることが確認できる．

補足：総和

Sympyでは総和も Sum を使って定義できる．この機能はフーリエ級数で少し登場する．

from sympy import Sum

i = symbols('i')
sum_expr = Sum(i, (i, 1, 10))
sum_expr

\[\displaystyle \sum_{i=1}^{10} i\]

総和の評価も doit() で計算できる．

sum_expr.doit()

\[\displaystyle 55\]