text_030_basic_calculation -- Sage

Hiroshi TAKEMOTO (take@pwv.co.jp)

数の計算・基本的な関数

よく使われる定数

数式でよく使われる定数をSageで使うには以下の様に表します。

hdr=["円周率", "自然体数の底", "虚数単位", "無限大"] sts=[pi, e, I, oo] html.table([(h, st, "%s"%st) for (h, st) in zip(hdr, sts)])





円周率
\pi
pi


自然体数の底
e
e


虚数単位
I
I


無限大
\infty
oo





円周率
\pi
pi


自然体数の底
e
e


虚数単位
I
I


無限大
\infty
oo

基本的な計算

四則演算をはじめ、Sageで使われる基本的な計算の表現方法を以下に示します。

# よく使われる表現 var('a b x n') hdr=["積", "商", "累乗", "平方根", "n乗根", "絶対値", "自然対数", "階乗"] sts=[a*b, a/b, a^b, sqrt(x), x^(1/n), abs(x), log(x), factorial(n)] html.table([(h, st, "%s"%st) for (h, st) in zip(hdr, sts)])





積
 $a b$ 
a*b


商
 $\frac{a}{b}$ 
a/b


累乗
 $a^{b}$ 
a^b


平方根
 $\sqrt{x}$ 
sqrt(x)


n乗根
 $x^{\left(\frac{1}{n}\right)}$ 
x^(1/n)


絶対値
 ${\left| x \right|}$ 
abs(x)


自然対数
 $\log\left(x\right)$ 
log(x)


階乗
 $n!$ 
factorial(n)





積

a*b


商

a/b


累乗

a^b


平方根

sqrt(x)


n乗根

x^(1/n)


絶対値

abs(x)


自然対数

log(x)


階乗

factorial(n)

1から10までの和は55ですが、これをリストとsum関数を使って計算すると、以下の様になります。

# 1から10までの和を計算 print 1+2+3+4+5+6+7+8+9+10 # リスト使って計算 L = range(1,11); print L print sum(L) # リストLの和を求める

55
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
55

55
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
55

リスト内包を使うと簡単にリストの要素を変更することができます。先ほどのリストLの要素を自乗に変えて、和を求めてみましょう。

	リスト内包の書式
	[ 式 for 変数 in リスト ]

L2 = [i^2 for i in L]; print L2 print sum(L2)

[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
385

[1, 4, 9, 16, 25, 36, 49, 64, 81, 100]
385

sumと並んでリストの要素の積を求めるprodもよく使われます。 1から5までの積をprodを使って計算してみましょう。

今回は、range関数の代わりにジェネレータを使って１から５までのリストを生成し、計算してみます。

	(1..5)

# 1から5までの積を計算 print 1*2*3*4*5 # リスト生成をまとめて print prod(1..5)

120
120

120
120

式の計算・簡単化

次に中学で習った式の因数分解や展開もSageを使うと簡単に計算することができます。

展開

以下の関数f1をSageを使って展開してみます。展開にはexpand関数を使います。 $$ f1(x) = (x - 1)(x^2 -1) $$

f1 = (x - 1)*(x^2 - 1); view(f1) # f1を定義 f2 = expand(f1); view(f2) # f1を展開し、f2に代入

 $\newcommand{\Bold}[1]{\mathbf{#1}}{\left(x^{2} - 1\right)} {\left(x - 1\right)}$ 
 $\newcommand{\Bold}[1]{\mathbf{#1}}x^{3} - x^{2} - x + 1$

多項式の係数は、coefficientsやcoeff関数で取得できます。

以下に、f2の係数をcoeeficents関数で取得する方法とxの2次の項の係数をcoeff関数で取得する方法を例として示します。

# 多項式の係数 print f2.coefficients(x) # 変数の[係数, 次数]のペアのリストを返す print f2.coeff(x, 2) # 変数と次数を指定して係数をを返す

[[1, 0], [-1, 1], [-1, 2], [1, 3]]
-1

[[1, 0], [-1, 1], [-1, 2], [1, 3]]
-1

因数分解

f1を展開した結果f2を因数分解してみましょう。因数分解にはfactor関数を使用します。

結果がf1と異なりますが、 $$ (x^2 - 1) = (x - 1)(x + 1) $$ の関係から、正しい結果になっていることが分かります。

view(factor(f2)) # f2を因数分解する

 $\newcommand{\Bold}[1]{\mathbf{#1}}{\left(x + 1\right)} {\left(x - 1\right)}^{2}$

部分分数分解

分数で表される式の積分やラブラス逆変換では部分分数分解を利用することがあります。

以下の様な関数f3を因数分解して、分数式f4に変えます。 $$ f3(x) = \frac{1}{(x - 2)} + \frac{1}{(x + 2)} $$

# 分数の因数分解 f3 = 1/(x+2)+1/(x-2); view(f3)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{x + 2} + \frac{1}{x - 2}$

式の部分分数分解には、prtial_fraction()メソッドを使います。f4にprtial_fractionメソッドを呼び出すとf3と同じ結果となります。

f4 = factor(f3); view(f4) # 分数式f3の因数分解 view(f4.partial_fraction()) # 部分分数分解

 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2 \, x}{{\left(x + 2\right)} {\left(x - 2\right)}}$ 
 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{x + 2} + \frac{1}{x - 2}$

式の整理

式を整理して簡単化する関数として、simplifyがあります。simplifyでは不要な項を消去するだけなので、 f1のような因数分解で整理されて式の場合には結果が変わりません。

さらに突っ込んだ整理をする場合には、simplify_fullメソッドを使います。

view(simplify(I + x -x)) # 簡単化で不要な項を消去 view(simplify(f1)) # 簡単化ではf1はそのまま view(f1.simplify_full()) # simplify_fullメソッドを使うと展開されて整理される

 $\newcommand{\Bold}[1]{\mathbf{#1}}i$ 
 $\newcommand{\Bold}[1]{\mathbf{#1}}{\left(x^{2} - 1\right)} {\left(x - 1\right)}$ 
 $\newcommand{\Bold}[1]{\mathbf{#1}}x^{3} - x^{2} - x + 1$

三角関数

三角関数も数式として計算するため、$ sin (\pi/4) $も数値ではなく、 $ \frac{1}{2} \sqrt{2} $の式が返ってきます。

値を得るには、N関数を使用します。

print sin(pi/4), N(sin(pi/4))

1/2*sqrt(2) 0.707106781186548

1/2*sqrt(2) 0.707106781186548

また三角関数への入力単位は度ではなく、ラジアンで指定します。

度とラジアンの変換には、以下の様な関数rad, degを定義すると便利です。

# 度とラジアンの変換関数 rad(x) = x*pi/180 deg(x) = x*180/pi

print sin(rad(45)), sin(pi/4) print atan(1), deg(atan(1))

1/2*sqrt(2) 1/2*sqrt(2)
1/4*pi 45

1/2*sqrt(2) 1/2*sqrt(2)
1/4*pi 45

三角関数の簡単化

三角関数や指数関数を含む式の簡単化には simplify_fullメソッドを使います。

例として、sin関数の倍角公式をsimplify_fullメソッドで求めてみます。

# 倍角公式 fs = sin(2*x); view(fs) view(fs.simplify_full())

 $\newcommand{\Bold}[1]{\mathbf{#1}}\sin\left(2 \, x\right)$ 
 $\newcommand{\Bold}[1]{\mathbf{#1}}2 \, \cos\left(x\right) \sin\left(x\right)$

更に、simplify_fullでは三角関数の公式を活用して式を整理します。

以下の例では、以下の三角関数の公式を使って簡単化しています。 $$ cos^2 x + sin^2 x = 1 $$

f4 = cos(x)^2-sin(x)^2; view(f4) # cos(x)^2 + sin(x)^2 = 1を使って view(f4.simplify_full()) # 簡単化する

 $\newcommand{\Bold}[1]{\mathbf{#1}}\cos\left(x\right)^{2} - \sin\left(x\right)^{2}$ 
 $\newcommand{\Bold}[1]{\mathbf{#1}}2 \, \cos\left(x\right)^{2} - 1$

導関数

高校で習った、導関数をSageを使って導いてみましょう。関数f(x)と平均変化率g(x)を以下の様に定義します。 $$ f(x) = \frac{1}{2} x^3 $$ $$ g(x)= \frac{f(x+h)−f(x)}{h} $$

変数x, hと関数f, gをSageで定義します。

x, h = var('x h') f(x) = x^3/2; view(f)

 $\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{1}{2} \, x^{3}$

# 平均変化率 g =(f(x + h) - f(x))/h; view(g)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(h + x\right)}^{3} - x^{3}}{2 \, h}$

gを展開して整理すると以下の様になり、h→0の極値（limit）をとった時、hが掛からないxの2次の項のみが残ります。

# 展開して整理すると g1 = g.simplify_full(); view(g1)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, h^{2} + \frac{3}{2} \, h x + \frac{3}{2} \, x^{2}$

求める導関数は、以下の様になります。

# h→0の極値を取ると導関数が求まる g2 = limit(g1, h=0); view(g2)

 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{3}{2} \, x^{2}$

特殊な関数のリミット

次のような特殊な関数fを考えます。 $$ f(x) = \frac{sin(x)}{x} $$

この関数にx=0を代入するとゼロ割のエラーとなりますが、limitを使ってx→0を求めると1となります。

Sageではグラフの表示でも極値をきちんと計算しているので、上記のf(x)もきちんと表示することができます。

# 特殊な関数のリミット f = sin(x)/x print f(x=0)

Traceback (click to the left of this block for traceback)
...
ValueError: power::eval(): division by zero

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_31.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("IyDnibnmrorjgarplqLmlbDjga7jg6rjg5/jg4Pjg4gKZiA9IHNpbih4KS94CnByaW50IGYoeD0wKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmph5BSYC/___code___.py", line 4, in <module>
    exec compile(u'print f(x=_sage_const_0 )
  File "", line 1, in <module>
    
  File "sage/symbolic/expression.pyx", line 4460, in sage.symbolic.expression.Expression.__call__ (build/cythonized/sage/symbolic/expression.cpp:25336)
  File "sage/symbolic/ring.pyx", line 788, in sage.symbolic.ring.SymbolicRing._call_element_ (build/cythonized/sage/symbolic/ring.cpp:9500)
  File "sage/symbolic/expression.pyx", line 4311, in sage.symbolic.expression.Expression.substitute (build/cythonized/sage/symbolic/expression.cpp:24455)
ValueError: power::eval(): division by zero

limit(f, x=0)

plot(f, [x, -10, 10])

関係式

Sageでの関係式の扱い方を以下の関係式eqを例に示します。 $$ x + 2 \le 0 $$

これをSageで表現すると以下の様になります。

eq = x + 2 <= 0

関係式から左辺、等号、不等号のオペレータ、左辺を取得するには、それぞれlhs, operator, rhs関数を使用します。

print eq.lhs() # 左辺 print eq.operator() # 等号、不等号などのオペレータ print eq.rhs() # 右辺

x + 2
<built-in function le>
0

x + 2
<built-in function le>
0

text_030_basic_calculation

3184 days ago by takepwave