2011/06/17からのアクセス回数 5356
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/7/
また、Sageのサーバを公開しているサイト(http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/)にユーザIDを作成することで、ダウンロードしたワークシートを アップロードし、実行したり、変更していろいろ動きを試すことができます。
高校で習う微分の公式は、以下の2つです。
sageで微分は、diff関数を使います。diff関数の使い方は以下の通りです。次数を省略すると1階微分となります。
diff(関数, 微分する変数, 次数)
上記公式をsageで実行すると以下のようになります。
sageへの入力:
# 変数と関数を生成 x, c, n = var('x c n') # 微分 print diff(c, x) print diff(x^n, x)
sageの出力:
0 n*x^(n - 1)
sageで変数xの関数fを以下のように定義します。
x = var('x') f = function('f', x)
この関数を微分すると、$f'$は以下のようにD[0](f)(x)のように返されます。
sageへの入力:
# 変数x, 関数fを定義 x = var('x') f = function('f', x) # 関数fを変数xで微分 diff(f, x)
sageの出力:
D[0](f)(x)
以下に示す微分の公式をSageで実行し、f'=D[0](f)(x)の置換を実行した結果を示します。
sageへの入力:
# df/dx = f', dg/dx = g'の置き換え関数の定義 repStr = lambda f : str(f).replace("D[0](f)", "f'").replace("D[0](g)", "g'")
#pre{{
}}
sageへの入力:
# 関数gを定義 g = function('g', x) # 微分公式をsageを使って実行 print repStr(diff(c*f, x)) print repStr(diff(f+g, x)) print repStr(diff(f*g,x)) print repStr(diff(f/g, x).simplify_full()) print repStr(diff(1/g,x)) print repStr(diff(f(g), x))
c*f'(x) f'(x) + g'(x) f(x)*g'(x) + g(x)*f'(x) -(f(x)*g'(x) - g(x)*f'(x))/g(x)^2 -g'(x)/g(x)^2 f'(g(x))*g'(x)
以下に主な関数を微分した結果を示します。
sageへの入力:
# いろんな関数の微分 var('x c n') sts=[c, x^n, sin(x), cos(x), exp(x), log(x)] html.table([("関数", "微分結果")] + [(st, diff(st, x)) for st in sts], header=true)
微分の応用として、接線の方程式とその法線を計算してみましょう。 関数f(x)を以下のように定義します。 $$ f(x) = x^3 - x^2 - 2*x $$ このとき、$x = x_1$での接線の式は、以下のようになります。 $$ y - y_{1} = f'(x_{1})(x - x_{1}) $$ また、法線は接線と直交することから、以下のようになります。 $$ y - y_{1} = -\frac{1}{f'(x_{1})}(x - x_{1}) $$ f(x)の微分した式をg(x)とし、x=0での接線と法線を計算し、プロットしてみます。
sageへの入力:
f(x) = x^3 - x^2 - 2*x; view(f)
sageへの入力:
g(x) = diff(f, x); view(g)
sageへの入力:
m = g(x=0); m
sageの出力:
-2
sageへの入力:
p1 = plot(f, [x, -2.5, 2.5]) p2 = plot(m*(x-0)+f(0), [x, -2.5, 2.5]) p3 = plot(-(x-0)/m+f(0),[x, -2.5, 2.5]) pt = point([0, f(0)]) (p1+p2+p3+pt).show(aspect_ratio=1, ymin=-2.5, ymax=2.5 )
sageでは積分は、integrate関数で計算します。
integrate(被積分関数, 積分変数) また integrate(被積分関数, 積分変数, 積分範囲) とすると定積分を計算します。
例として、\(\int x dx\)を計算した後、その微分を求めています。 元のxが返されています。
sageへの入力:
f = integrate(x, x); view(f)
sageへの入力:
diff(f)
sageの出力:
x
以下に主な関数を積分した結果を示します。
sageへの入力:
# いろんな関数の積分 var('x c') f = function('f', x) sts=[c, 1/x, exp(x), log(x), sin(x), diff(f, x)/f(x)] html.table([("関数", "積分結果")] + [(repStr(st), integrate(st, x)) for st in sts], header=true)
次に定積分$\int_{0}^{\pi/2}sin(x) dx$の例を示します。
$sin(x)$を0から$\pi/2$の範囲でプロットすると以下のようになります。
sageへの入力:
plot(sin, [x, 0, pi/2], fill=True)
sinの積分は、-cosですから、定積分は以下のようになります。 $$ [ -cos(x) ]_0^\pi = \{ 0 - (-1) \} = 1 $$
それでは、Sageを使って、\(sin(x)\)を0から\(\pi/2\)の範囲で積分してみましょう。
sageへの入力:
integrate(sin(x), x, 0, pi/2)
sageの出力:
1
もう一つ、以下の関数を積分した例をみてみましょう。 $$ \int_{0}^{3}x^2 + 2sin(x) dx $$
先ほどと同様に、積分範囲をプロットしてみます。
sageへの入力:
plot(x^2 + 2*sin(x), [x, 0, 3], fill=True)
積分の結果は、関数で返されるのでN(_)を使って数値に変換します。
このように関数のプロット結果から推定される概算と積分の結果をみながら、 計算結果を確認しながら、進められるのもSageの特徴です。
sageへの入力:
integrate(x^2+2*sin(x), x, 0, 3)
sageの出力:
-2*cos(3) + 11
sageへの入力:
N(_)
sageの出力:
12.9799849932009
すべての関数が定積分できるとは限りません、そこでSageには数値積分を行う numerical_integral関数が用意されています。
numerical_integral関数の使い方は、以下の通りです。 関数の戻り値は、積分結果と誤差のタプルが返されます。
numerical_integral(被関分館数, 積分範囲)
例として、\(sin(x)\)を0から\(\pi/2\)で数値積分した結果を以下に示します。
sageへの入力:
sigma, error = numerical_integral(sin(x), 0, pi/2); (sigma, error)
(1.0, 1.1102230246251565e-14)
積分の応用例として、サイクロイドの計算をしてみます。
サイクロイドは、tをパラメータとして、以下のような式で表されます。 $$ \left\{\begin{eqnarray} x &=& 2(t-sin(t)) \\ y &=& 2(1-cos(t)) \end{eqnarray}\right. $$
この曲線の長さは、以下の積分で計算することができます。 $$ \int_{0}^{2\pi}\sqrt{(\frac{dx}{dt})^2+(\frac{dy}{dt})^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))
sageへの入力:
cycloid = sqrt(diff(x,t)^2+diff(y,t)^2)
残念ながら、integrate関数を使うと符号がマイナスの原因が不明ですが、-16が返ってきます。
sageへの入力:
integrate(cycloid, t, 0, 2*pi)
sageの出力:
-16
次に数値積分で計算してみます。16に近い値が求まります。
Sageのような数式処理で計算する場合には、途中計算のチェックや他の手法との相互 チェックをするように心がけるとよいでしょう。
sageへの入力:
numerical_integral(cycloid, 0, 2*pi)
sageの出力:
(15.999999999999998, 1.7763568394002502e-13)
皆様のご意見、ご希望をお待ちしております。