1048-計算してみよう.sagews

Hiroshi TAKEMOTO([email protected])

計算してみよう

変数の宣言と値の代入

これまで式の変数は宣言しないで使っていましたが、正式には変数はvar関数で宣言して 使います。 変数の値を指定して、関数の返す値を調べるには、

f1(x=1.5)

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

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

0.625000000000000

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

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

7

規則の代入

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

f6=x^2 + 1; f6

x^2 + 1

f6.subs_expr(x^2 == w)

w + 1

関数の定義

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

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

です。

python関数の定義

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

pythonでは

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

のように定義します。

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

def cube(x):
    return x^3

cube(2)

8

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

(x + y)^3

lambda関数の定義

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

lambda 変数のリスト: 式

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

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$を代入して解の値を取得します。

# 方程式の解法
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$をセットするようにします。こうすれば簡単に解の値が確認できます。

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$の範囲で数値解析します。

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

-0.5

多項式の解

一次方程式ではあまりに簡単なので、多項式の解を求めてみます。

$f(x)=x^3-2x+4$

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

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

plot(f, -2.5, 2.5)

solveの結果、$x=-2, x=1-i, x=1+i$を得ました。

solve(f, x)

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

因数分解の結果からも、この解が正しいことがわかります。

factor(f)

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

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

find_rootの値は、以下の通りです。

find_root(f, -3, 3)

-1.9999999999999949

連立方程式の解

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

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

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

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]]

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

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{darray}{rcl} x^2 +y^2 &=& 1 \\ y &=& x -1 \end{darray}$$

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

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としました。

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) = x^3$とします。

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

# 導関数
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)$を展開し、整理すると

g1= g.rational_simplify(); g1

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

となります。

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

limit(g1, h=0)

3/2*x^2

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

diff(f(x), x)

3/2*x^2

微分

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

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

です。

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

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

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

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

diff(c, x)

0

diff(x^n, x)

n*x^(n - 1)

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

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

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

diff(f, x)

D[0](f)(x)

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

• $(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)

となり、微分の各公式を表しています。

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)

いろんな関数の微分

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

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

# いろんな関数の微分
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'(x_{1})}(x - x_{1})$

となります。

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

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 x dx$を計算した後、その微分を求めています。

f = integrate(x); f

1/2*x^2

diff(f)

x

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

• $\int c dx$ 定数$c$
• $\int 1/x dx$
• $\int e^x dx$
• $\int log(x) dx$
• $\int sin(x) dx$
• $\int f'(x)/f(x) dx$

# いろんな関数の積分
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$の計算を示します。

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

1

もう一つ

$\int_{0}^{3}x^2 + 2sin(x) dx$

の例を示し、

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

-2*cos(3) + 11

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

N(_)

12.9799849932009

数値積分

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

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

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

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

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

(1.9999999999999998, 2.2204460492503128e-14)

[removed]

積分の応用

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

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

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

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

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

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

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)