数式の変形は,手で直すほうが圧倒的に早くきれいになる場合が多い.しかし,テイラー展開や,複雑な積分公式,三角関数とexp関数の変換などの手間がかかるところを,数式処理ソフトは間違いなく変形してくれる.ここで示すコマンドを全て覚える必要は全くない.というか忘れるもの.ここでは,できるだけコンパクトにまとめて,悩んだときに参照できるようにする.初めての人は,ざっと眺めた後,鉄則からじっくりフォローせよ.
from sympy import *
x,y=symbols('x y')
simplify(3*x+4*x+2*y)
7*x + 2*y
exp1=3*sin(x)**3-sin(x)*cos(x)**2
simplify(exp1)
4*sin(x)**3 - sin(x)
expand((x+y)**2)
x**2 + 2*x*y + y**2
factor(4*x**2-6*x*y+2*y**2)
2*(x - y)*(2*x - y)
simplify((x+y)/(x**2-3*x*y-4*y**2))
1/(x - 4*y)
together(1/x+1/y)
(x + y)/(x*y)
?sqrtdenest
a,b,c = symbols('a b c')
collect(4*a*x**2-3*y**2/x+6*b*x*y+3*c*y+2*y**2,y)
4*a*x**2 + y**2*(2 - 3/x) + y*(6*b*x + 3*c)
分数を英語でfractionというが,これでとり出せる.出力は,numerator(分子),denominator(分母)の順.
a,b,c, x = symbols('a b c x')
frac = (1+x)**2/(a*(x-1)+b*(x+3))
n,d = fraction(frac)
pprint(n)
pprint(d)
2 (x + 1) a⋅(x - 1) + b⋅(x + 3)
eq = a+x+b*x+c*x**2
pprint(eq)
eq.coeff(x)
2 a + b⋅x + c⋅x + x
b + 1
eq.as_coefficients_dict()
defaultdict(int, {a: 1, x: 1, b*x: 1, c*x**2: 1})
eq0 = Poly(eq, x)
coeffs = eq0.coeffs()
print(coeffs[::-1])
[a, b + 1, c]
方程式などで,係数として分数を使うとうまく扱ってくれません. そんなときは,Rationalで有理数と指定してしまいます.
factor(1/4*x**2-y**2)
1.0*(0.25*x - 0.5*y)*(0.25*x + 0.5*y)
factor(Rational(1,4)*x**2-y**2)
(x - 2*y)*(x + 2*y)/4
factor(x**3-8)
c = symbols('c')
factor(a*b*c - a*b*x - a*c*x + a*x**2 + b*c*x - b*x**2 - c*x**2 + x**3)
apart((x+1)/((x-1)*(x**2+1)**2))
apart(3/(1-x**4))
from sympy import *
radsimp(8/(3-sqrt(5)))
radsimp(2/(2+sqrt(5)))
simplify(8/(3-sqrt(5))-2/(2+sqrt(5)))
10
k,x = symbols('k x')
eq=x**2 +2*k*x+5-k
sol = solve(eq,x)
solve(sol[0]-sol[1],k)
[-1/2 + sqrt(21)/2, -sqrt(21)/2 - 1/2]
どうしても解かなければならない課題を前にコマンドリファレンスのあちこちを参照しながら解いていくのが数式処理を修得する最速法である.とびかかる前にちょっとした共通 のコツがある.それをここでは示す.数式処理ソフトでの数式処理とは,数式処理ソフトが『自動的にやって』くれるのではなく,実際に紙と鉛筆で解いていく手順を数式処理ソ フトに『やらせる』ことであることを肝に銘じよ.
数式処理ソフトの習得にあたって初心者がつまづく共通の過ちを回避する鉄則がある.
# kernel->Restart
from sympy import *
init_session()
として,sympyのimport,init_sessionによる綺麗な表示を強制しておく.
ValueError: The same variable should be used in
all univariate expressions being plotted.
今まで出てきたコマンドを使えば,典型的なセンター試験の問題を解くのも容易である.以下の例題を参照して課題を解いてみよ.使うコマンドは,unapply, solve, diff, expand(展開), factor(因数分解)とsubs(一時的代入)である.
$a,b$を定数とし, $a \neq 0$とする.2次関数
$$ y = a\,x^2-b\,x-a+b\,\cdots(1) $$のグラフが点 $(-2, 6)$ を通るとする.このとき
$$ b = -a+\fbox{ ア } $$であり,グラフの頂点の座標を$a$を用いて表すと
$$ \left(\frac{-a+\fbox{ イ }}{\fbox{ ウ }\,a}, -\frac{(\fbox{ エ }\,a- \fbox{ オ })^2}{\fbox{ カ }\,a}\right) $$である (2008 年度大学入試センター試験数学 I 第2問より抜粋).
from sympy import *
init_session()
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b
bb = solve(f.subs({x:-2})-6,b)[0]
pprint(bb)
x0=solve(diff(f,x),x)[0].subs({b:bb})
pprint(x0)
y0=factor(f.subs({x:x0,b:bb}))
pprint(y0)
IPython console for SymPy 1.0 (Python 3.6.1-64-bit) (ground types: python) These commands were executed: >>> from __future__ import division >>> from sympy import * >>> x, y, z, t = symbols('x y z t') >>> k, m, n = symbols('k m n', integer=True) >>> f, g, h = symbols('f g h', cls=Function) >>> init_printing() Documentation can be found at http://docs.sympy.org/1.0/ -a + 2 -a + 2 ────── 2⋅a 2 -(3⋅a - 2) ──────────── 4⋅a
解答の一例は以上のとおりであるが,これはだいぶ整理されたあとの記述になる. 普通の答案ではここまで綺麗に初めから求まるものではない.
そこで,綺麗な解答を導く一つ一つのステップをはしょることなくここで示しておく. 今後巡り会うであろう解答例やコードサンプルにおいても, ここで示した導出のコンセプトを参考に数式処理コードを理解していってもらいたい.
まずはRestartをかける.これもcodeに示しておくと忘れない. それとsympyのimportとinit_sessionによる綺麗なprint outを癖付けておく.
# kernel->Restart
from sympy import *
init_session()
IPython console for SymPy 1.0 (Python 3.6.1-64-bit) (ground types: python) These commands were executed: >>> from __future__ import division >>> from sympy import * >>> x, y, z, t = symbols('x y z t') >>> k, m, n = symbols('k m n', integer=True) >>> f, g, h = symbols('f g h', cls=Function) >>> init_printing() Documentation can be found at http://docs.sympy.org/1.0/
与えられた2次関数を$f(x)$で定義する
f = a*x**2-b*x-a+b
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-2-de84ddc8e8d5> in <module>() ----> 1 f = a*x**2-b*x-a+b NameError: name 'a' is not defined
すると早速エラー.sympyでは変数を明示的に宣言しないといけない.
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b # 与えられた2次関数を$f(x)$で定義する
関数fが点(−2,6)を通る式を立てる.まずはx=-2を代入(subs)
f.subs({x:-2})
こいつがy=-6をみたすのだから,yに前の式を代入して,左辺にまとめて見る.
y = f.subs({x:-2})
y-6
これをeq1と置いて,bについて解く.
eq1 = y-6
solve(eq1,b)
これがbの表式.配列の一番目の要素となっているので,それを取り出しておく.
bb = solve(eq1,b)[0]
pprint(bb)
-a + 2
頂点の座標(x0,y0)で傾きが0になることを用いて解いていく.鉄則3を思い出して,中から見ていく. 微分diffを思い出して,
diff(f,x)
これをxについて解く.
solve(diff(f,x),x)
solve(diff(f,x),x).subs({b:bb})
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) <ipython-input-10-33528e910811> in <module>() ----> 1 solve(diff(f,x),x).subs({b:bb}) AttributeError: 'list' object has no attribute 'subs'
oops. このエラーはlistに対してはsubsは使えないとのこと. 先に中身を取り出しておいて,代入.
solve(diff(f,x),x)[0].subs({b:bb})
これが頂点x0.これをfに代入してやる.ついでにb=bbも...
f.subs({x:x0,b:bb})
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-12-0738cd0e61b9> in <module>() ----> 1 f.subs({x:x0,b:bb}) NameError: name 'x0' is not defined
せやった.x0に代入してない.
x0=solve(diff(f,x),x)[0].subs({b:bb})
f.subs({x:x0,b:bb})
こいつを簡単化simplifyすればOK.
simplify(f.subs({x:x0,b:bb}))
あかん.ここは明示的に因数分解factorを指定する.
factor(f.subs({x:x0,b:bb}))
でや?! これをまとめて表示すると,
# kernel->Restart
from sympy import *
init_session()
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b # 与えられた2次関数を$f(x)$で定義する
y = f.subs({x:-2})
eq1 = y-6
bb = solve(eq1,b)[0]
pprint(bb)
x0=solve(diff(f,x),x)[0].subs({b:bb})
pprint(x0)
factor(f.subs({x:x0,b:bb}))
IPython console for SymPy 1.0 (Python 3.6.1-64-bit) (ground types: python) These commands were executed: >>> from __future__ import division >>> from sympy import * >>> x, y, z, t = symbols('x y z t') >>> k, m, n = symbols('k m n', integer=True) >>> f, g, h = symbols('f g h', cls=Function) >>> init_printing() Documentation can be found at http://docs.sympy.org/1.0/ -a + 2 -a + 2 ────── 2⋅a
あれれ...なんか途中が抜けとるで??
途中の重要な結果は変数に置き換えて明示pprintしておくのが鉄則1. それを意識して修正すると,
# kernel->Restart
from sympy import *
init_session()
a,b,x = symbols('a b x')
f = a*x**2-b*x-a+b # 与えられた2次関数を$f(x)$で定義する
bb = solve(f.subs({x:-2})-6,b)[0]
pprint(bb)
x0=solve(diff(f,x),x)[0].subs({b:bb})
pprint(x0)
y0=factor(f.subs({x:x0,b:bb}))
pprint(y0)
IPython console for SymPy 1.0 (Python 3.6.1-64-bit) (ground types: python) These commands were executed: >>> from __future__ import division >>> from sympy import * >>> x, y, z, t = symbols('x y z t') >>> k, m, n = symbols('k m n', integer=True) >>> f, g, h = symbols('f g h', cls=Function) >>> init_printing() Documentation can be found at http://docs.sympy.org/1.0/ -a + 2 -a + 2 ────── 2⋅a 2 -(3⋅a - 2) ──────────── 4⋅a
(例題1に続いて)
さらに,2次関数(1)のグラフの頂点の$y$座標が-2であるとする.このとき,$a$は \begin{equation*} \fbox{ キ }\,a^2-\fbox{ クケ }\,a+\fbox{ コ } = 0 \end{equation*} を満たす.これより,$a$の値は \begin{equation*} a = \fbox{ サ }, \frac{\fbox{ シ }}{\fbox{ ス }} \end{equation*} である. 以下,$a = \frac{\fbox{ シ }}{\fbox{ ス }}$であるとする.
このとき,2次関数(1)のグラフの頂点のx座標は$\fbox{ セ }$であり,(1)のグラフと$x$軸の2交点の$x$座標は$\fbox{ ソ },\fbox{ タ }$である.
また,関数(1)は$0 \leqq x \leqq 9$において
$x = \fbox{ チ }$のとき,最小値$\fbox{ ツテ }$をとり,
$x = \fbox{ ト }$のとき,最大値$\frac{\fbox{ ナニ }}{\fbox{ ヌ }}$をとる.
(2008 年度大学入試センター試験数学 I より抜粋).
# 実践の解答例に続けて...
eq2 = expand((y0 +2)*(-4*a)) # -4aをかけて,展開
pprint(eq2) # キクケコ
s_a=solve(eq2,a) # 方程式を解く
pprint(s_a) # サシス
aa = min(s_a) #二つの解の小さい方が題意,シス
pprint(aa)
pprint(x0.subs({a:aa})) # セ
eq3 = f.subs({b:bb}).subs({a:aa}) # fに得られたa,bを代入
pprint(solve(eq3,x)) # そいつを解いてx軸との交点を,ソタ
2 9⋅a - 20⋅a + 4 [2/9, 2] 2/9 4 [1, 7]
%matplotlib inline
plot(eq3, (x,0,9)) # 得られたeq3をplot
<sympy.plotting.plot.Plot at 0x110133f98>
print(eq3.subs({x:4})) # min at x=4
print(eq3.subs({x:9})) # 目視でmax at x=9
-2 32/9
$P = x(x+3)(2x-3)$とする. また,$a$を定数とする. $x = a+1$のときの $P$の値は \begin{equation*} 2a^3+\fbox{ ア }a^2+\fbox{ イ }a-\fbox{ ウ } \end{equation*} である.
$x=a+1$のときの$P$の値と,$x=a$のときの$P$の値が等しいとする.このとき,$a$は \begin{equation*} 3a^2+\fbox{ エ }a-\fbox{ オ } = 0 \end{equation*} を満たす.したがって \begin{equation*} a = \frac{\fbox{ カキ }\pm \sqrt{\fbox{ クケ }}}{\fbox{ コ }} \end{equation*} である.
(2008 年度大学入試センター試験数学 I 第4問より抜粋)
from sympy import *
a,b,x = symbols('a b x')
P = x*(x+3)*(2*x-3) # $P$を$x$の関数として定義
pprint(expand(P.subs({x:a+1}))) # $$P(a)$を形式的に出してみる.
eq1 = expand(P.subs({x:a+1}) - P.subs({x:a})) #2式を差し引く
pprint(eq1/2) # 出題にそろえるため2で割っている.
solve(eq1,a) # eq1=0をaについて解く(solve)
3 2 2⋅a + 9⋅a + 3⋅a - 4 2 3⋅a + 6⋅a - 2
数式処理ソフトで実際に数式をいじる状況というのは,ほとんどの場合が既知の数式変形のフォローだろう.例えば,論文で「(1)式から(2)式への変形は自明である」とかいう 文章で済ましている変形が本当にあっているのかを確かめたい時.一番単純なやり方は自明と言われた前後の式が一致していることを確かめるだけで十分である. 最も単純な確認法は以下の通り,変形の前後の式を手入力してその差をexpandした結果が0か否かでする.
from sympy import *
x = symbols('x')
ex1=(x-3)**4
ex2 = x**4-12*x**3+54*x**2-108*x+81
expand(ex1-ex2)
0ならば式の変形は保証されているので,その導出が間違いでなく誤植などもないことが確認できる.ただ,これだけでは変形の哲学や技法が身に付くわけではない.あくまでも 苦し紛れのデフォルトであることは心に留めておくように. 論理値(True or False)として確かめたいときには,==0で出てくる.
以下に示す積分を実行せよ.
$$ \int _{-\infty }^{\infty }x \exp(-\beta c\,{x}^{2}) \left( 1+\beta g\,{x}^{3} \right) {dx} $$最新版のMapleでは改良が施されていて,このような複雑な積分も一発で求まる.
maple
> int(f1(x),x=-infinity..infinity);
$$
\left\{\, \begin {array}{cc} { \frac{3}{4}\,\frac {g \sqrt{\pi }}{\beta\,{c}^{2} \sqrt{c\beta}}}&csgn \left( c\beta \right) =1\\ \infty&otherwise\end {array} \right.
$$ここでは,$c \beta$が正の場合(csgn(beta c)=1)とそれ以外の場合(otherwise)に分けて答えを返している.しかしこのような意図したきれいな結果をいつも数式処理ソフトが返してくれるとは限らない.これだけしか知らないと,なにかうまくいかないときにお手上げになってしまう.
from sympy import *
init_printing()
x,beta,c,g = symbols('x beta c g')
f1 = x*exp(-beta*c*x**2)*(1+beta*g*x**3)
pprint(f1)
pprint(integrate(f1,(x,-oo,oo)))
2 ⎛ 3 ⎞ -β⋅c⋅x x⋅⎝β⋅g⋅x + 1⎠⋅ℯ ⎧ 3⋅√π⋅g ⎪───────────────────────────────────── for │periodic_argument(polar_lift(β)⋅p ⎪ 3/2 3/2 ⎪4⋅c⋅polar_lift (β)⋅polar_lift (c) ⎪ ⎪ ∞ ⎨ ⌠ ⎪ ⎮ 2 ⎪ ⎮ ⎛ 3 ⎞ -β⋅c⋅x ⎪ ⎮ x⋅⎝β⋅g⋅x + 1⎠⋅ℯ dx otherwise ⎪ ⌡ ⎪ -∞ ⎩ π olar_lift(c), ∞)│ < ─ 2
この辺りの記述をやり直したけど, assumeのためにassumpitons moduleを使ってやったんだけど, どうもうまくいかない.
assumes = Q.positive(beta), Q.positive(c)
with assuming(*assumes):
print ask(beta*c, positive)
とかは通るけど,integrateの結果は変わらない. 結局,symbols('x',positive = True)とかで初めに指定するとすっと通る.