sage/計算してみよう

FrontPage

2009/11/14からのアクセス回数 12281

このページで紹介したsageノートブックは、以下のURLから参照することができます。

http://www.sagenb.org/home/pub/1048/

変数の宣言と値の代入 †

sageの式で変数として認識させるには変数をvar関数で宣言しなくてなりません。

変数の宣言は、

var('変数名')

複数の変数を宣言する場合には、スペースを空けて指定します。 宣言される変数を参照したりする場合には、

x, y = var('x y')

のようにpython変数x, yに宣言したsage変数を代入します。

変数の値を指定して、関数の返す値を調べるには、

f1(x=1.5)

のように式を表す変数にカッコを付けて、x=1.5と値を指定するとその時の関数の値が表示されます。

sage入力：

x = var('x')
f1 = (x - 1)*(x^2 - 1)
f1(x=1.5)

0.625000000000000

変数をx, yと複数存在する場合の例をです。f5に\(f5(x,y)=x2+y\)を定義し、f5(x=2, y=3)でx=2, y=3の時の値を出力します。

sage入力：

var('x y w')
f5=x^2 + y
f5(x=2, y=3)

7

規則の代入 †

規則の代入には、subs_exprを使います。

sage入力：

f6=x^2 + 1; f6

x^2 + 1

sage入力：

f6.subs_expr(x^2 == w)

w + 1

関数の定義 †

sageで関数を定義する方法が２つあります。

• pythonの関数を定義する
• lambda式を使って関数を定義する

です。

python関数の定義 †

手続き的な関数を定義する場合には、python関数を使うと便利です。

pythonでは

def 関数名(引数のリスト):
    関数の本体

のように定義します。

以下にx3を返す関数cubeをpythonの関数を使って定義しました。

sage入力：

def cube(x):
    return x^3

cube(2)

8

x, y = var('x y')
cube(x+y)

(x + y)^3

lambda関数の定義 †

つぎにlambda式を使って関数を定義する方法を示します。 x3のように式で表現できる場合には、lambda関数が便利です。

lambda 変数のリスト: 式

のように定義します。 先ほどと同じx3を返す関数の例を以下に示します。

sage入力：

f = lambda x: x^3
f(2)

8

f(x+y)

(x + y)^3

方程式の解法 †

数式処理で助かる機能は、方程式の解法でしょう。

方程式の解法は、solve関数を使って以下のように行います。

solve(解放したい関数また関数のリスト, 解を得る変数)

例として、\(f(x)=ax+b\)の一次方程式を解いてみます。結果から\(x=−b/a\) が得られたので、\(a=2,b=1\)を代入して解の値を取得します。

sage入力：

# 方程式の解法
var('x a b')
f = a*x + b

sln = solve(f, x); sln

[x == -b/a]

x = -b/a
x(a=2, b=1)

-1/2

これだと、解から式をもう一度入力しなくてはなりません。そこで、a,bも変数として関数を定義し、solveでa=2,b=1をセットするようにします。こうすれば簡単に解の値が確認できます。

sage入力：

var('x a b')
g = lambda x, a, b : a*x + b

solve(g(x, a, b), x)

[x == -b/a]

solve(g(x, a=2, b=1), x)

[x == (-1/2)]

解の値がほしいときには、find_root関数で数値解析で値を計算します。

find_root(関数, 計算する変数の開始値、終了値)

先ほどの関数g(x):a=2,b=1の解をx=−1,1 の範囲で数値解析します。

sage入力：

find_root(g(x, a=2, b=1), -1, 1)

-0.5

多項式の解 †

一次方程式ではあまりに簡単なので、多項式の解を求めてみます。 $$ f(x)=x^3−2x+4 $$

のプロット図を以下に表示します。x軸と交差しているのはx=−2 です。-1と1付近に極値があります。

sage入力：

var('x')
f = x^3-2*x+4
plot(f, -2.5, 2.5)

solveの結果、\(x=−2,x=1−i,x=1+i \)を得ました。 sage入力：

solve(f, x)

[x == (-I + 1), x == (I + 1), x == -2]

sage入力：

factor(f)

(x + 2)*(x^2 - 2*x + 2)

plotで可視化することでfind_rootでの計算範囲を大まかに把握することができます。

図から-3から3の間に解があることが見て取れたので、この範囲で数値解を求めます。 sage入力：

find_root(f, -3, 3)

-1.9999999999999949

連立方程式の解 †

同様にsolveに関数に式のリストを与えることで、連立方程式の解を得ることができます。

$$ \begin{eqnarray} x^2 - 2y &=& 2 \\ x - y &=& 1 \end{eqnarray} $$

の解を例に以下に示します。

sage入力：

var('x y')
f = [y == 1/2*x^2 - 1, y == x - 1]
solve(f, x, y)

[ [x == 2, y == 1], [x == 0, y == -1] ]

上記関数のプロット結果は、以下のとおりです。

sage入力：

p1 = plot(lambda x : 1/2*x^2 - 1, (x, -1, 3))
p2 = plot(lambda x : x - 1, (x, -1, 3))
(p1+p2).show()

もう一つ、

$$ \begin{eqnarray} x^2 + y^2 &=& 1 \\ y &=& x - 1 \end{eqnarray} $$

の解を例に以下に示します。

sage入力：

var('x y')
f = [x^2 + y^2 == 1, y == x - 1]
solve(f, x, y)

[ [x == 1, y == 0], [x == 0, y == -1] ]

この関数をimplicit_plotを使って関係式の形をそのまま使ってプロットしてみます。 円のなので、aspect_ratio=1としました。

sage入力：

var('x y')
p1 = implicit_plot(x^2 + y^2 == 1, (x, -1, 1), (y, -1, 1))
p2 = plot(lambda x : x - 1, (x, -1, 1))
(p1+p2).show(aspect_ratio=1)

導関数 †

高校で習った、導関数をsageを使って導いてみましょう。

関数\(f(x)=\frac{1}{2} x^3\)とします。

平均変化率は、\(g(x)=hf(x+h)−f(x)\)を以下のように定義します。

sage入力：

# 導関数
def f(x):
    return (x^3)/2

var('x, h')
g=(f(x + h) - f(x))/h; g

1/2*((h + x)^3 - x^3)/h

g(x)を展開し、整理すると

sage入力：

g1= g.rational_simplify(); g1

1/2*h^2 + 3/2*h*x + 3/2*x^2

となります。

ここで、h→0 の極値を取ったものが求める導関数\(f′(x)\)です。

sage入力：

limit(g1, h=0)

3/2*x^2

結果は、f(x)をxで微分したものと同じになります。

sage入力：

diff(f(x), x)

3/2*x^2

微分 †

高校で習う微分の公式は、

• \(\frac{d}{dx} c = 0\)
• \(\frac{d}{dx} x^n = n x^{x-1}\)

です。

sageの微分diff関数は、以下のように使います。

diff(関数, 微分する変数)

上記公式のsageで結果は、以下のようになります。

sage入力：

# 微分
x, c, n = var('x c n')
f = function('f', x)
g = function('g', x)

print diff(c, x)
print diff(x^n, x)

0

n*x^(n - 1)

sageで変数xの関数fを以下のように定義し微分すると、\(f′\)は以下のようにD[0](f)(x)となります。

x = var('x')
f = function('f', x)

これで変数xの関数fを定義したことになります。

sage入力：

diff(f, x)

D[0](f)(x)

このD[0](f)を\(f′\)と読み替えると、

sage入力：

print diff(c*f, x)
print diff(f+g, x)
print diff(f*g,x)
print diff(f/g, x).simplify_full()
print diff(1/g,x)
print diff(f(g), x)

c*D[0](f)(x)

D[0](f)(x) + D[0](g)(x)

D[0](f)(x)*g(x) + f(x)*D[0](g)(x)

(D[0](f)(x)*g(x) - D[0](g)(x)*f(x))/g(x)^2

-D[0](g)(x)/g(x)^2

D[0](f)(g(x))*D[0](g)(x)

は、以下のようようになり、微分の各公式を表しています。

• \( (cf(x))′=cf′(x)\)
• \( (f(x)+g(x))′=f′(x)+g′(x)\)
• \( (f(x)g(x))′=f′(x)g(x)+f(x)g′(x)\)
• \( f(x)/g(x)=(f(x)′g(x)−g′(x)f(x))/g(x)^2 \)
• \( (1/g(x))′=−g′(x)/g(x)^2\)
• \( (f(g(x)))′=f′(g(x))g′(x)\)

いろんな関数の微分 †

以下の関数を微分した結果を示します。

• \( (x^2+1)^3\)
• \( sin(x)^3\)
• \( sin^{-1}(x)\)
• \( e^x \)
• \( log(x)\)

sage入力：

# いろんな関数の微分
print diff((x^2+1)^3).factor()
print diff(sin(x)^3, x)
print diff(arcsin(x), x)
print diff(exp(x), x)
print diff(log(x), x)

6*(x^2 + 1)^2*x

3*sin(x)^2*cos(x)

1/sqrt(-x^2 + 1)

e^x

1/x

微分の応用 †

微分の応用として、接線の方程式とその法線を計算してみましょう。\( f(x)=e^{−x^2}\)の接線は、

$$ y−y_1=f′(x_1)(x−x_1) $$

となります。 また、法線は、 $$ y−y_1=\frac{1}{f′(x1)}(x−x_1) $$ となります。

微分した式をgとし、x=1での接線と法線を計算し、プロットしてみます。

sage入力：

f = exp(-(x^2)); f

e^(-x^2)

g = diff(f, x); g

-2*x*e^(-x^2)

m = g(x=1); m

-2*e^(-1)

p1 = plot(f, (x, -2, 2))
p2 = plot(m*(x-1)+f(1), (x, -2, 2))
p3 = plot(-(x-1)/m+f(1), (x, -2, 2))
pt = point([1, f(1)])
(p1+p2+p3+pt).show(aspect_ratio=1, ymin=-1, ymax=2)

積分 †

sageでは積分は、integrate関数で計算します。

integrate(被積分関数, 積分変数)
また
integrate(被積分関数, 積分変数, 積分範囲)

とすると定積分を計算します。

\(\int xdx\)を計算した後、その微分を求めています。

sage入力：

f = integrate(x); f

1/2*x^2

diff(f)

1

以下の関数を積分した結果示します。

• \( \int cdx \)定数c
• \( \int \frac{1}{x}dx\)
• \( \int e^x dx\)
• \( \int log(x)dx \)
• \( \int sin(x)dx \)
• \( \int \frac{f'(x)}{f(x)}dx\)

sage入力：

# いろんな関数の積分
x = var('x')
c = var('c')
f = function('f', x)
g = function('g', x)

print integrate(c, x)
print integrate(1/x, x)
print integrate(exp(x), x)
print integrate(log(x), x)
print integrate(sin(x), x)
print integrate(diff(f, x)/f(x), x)

c*x

log(x)

e^x

x*log(x) - x

-cos(x)

log(f(x))

定積分 †

次に定積分\( \int_0^{\pi/2}sin(x)dx \) の計算を示します。

sage入力：

integrate(sin(x), x, 0, pi/2)

1

もう一つ $$ \int_0^3 x^2 + 2 sind(x) dx $$

の例を示し、

sage入力：

integrate(x^2+2*sin(x), x, 0, 3)

-2*cos(3) + 11

上記結果を数値とすると以下のようになります。

N(_)

12.9799849932009

数値積分 †

数値積分は、numerical_integral関数を使って計算します。

numerical_integral(被関分館数, 積分範囲)

で戻り値は、積分結果と誤差を返します。

以下は、\(sin(x\))を0から\(\pi\)で数値積分した結果です。

sage入力：

sigma, error = numerical_integral(sin(x),  0, pi); (sigma, error)

(1.9999999999999998, 2.2204460492503128e-14)

積分の応用 †

積分の応用例として、サイクロイドの計算をしてみます。

サイクロイドは、tをパラメータとして、以下のような式で表されます。

$$ \begin{eqnarray} x &=& 2(t-sin(t)) \\ y &=& 2(1-cos(t)) \end{eqnarray} $$

この曲線の長さを計算すると、

$$ \int_{0}^{2\pi}\sqrt{\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2}dt $$

となります。解析解では、この結果は16となります。

sage入力：

t = var('t')
x = 2*(t-sin(t))
y = 2*(1-cos(t))
parametric_plot([x, y], (t, 0, 2*pi))

cycloid = sqrt(diff(x,t)^2+diff(y,t)^2)

残念ながら、定積分は結果がでませんでした。

integrate(cycloid, t, 0, 2*pi)

integrate(sqrt(4*(cos(t) - 1)^2 + 4*sin(t)^2), t, 0, 2*pi)

数値積分で、期待した結果となりました。

numerical_integral(cycloid, 0, 2*pi)

(15.999999999999998, 1.7763568394002502e-13)

コメント †

選択肢	投票
おもしろかった	16
そうでもない	0
わかりずらい	1

皆様のご意見、ご希望をお待ちしております。

お名前:

(Input image string)