#freeze [[FrontPage]] #contents 2009/11/14からのアクセス回数 &counter; このページで紹介したsageノートブックは、以下のURLから参照することができます。 http://www.sagenb.org/home/pub/1048/ ** 変数の宣言と値の代入 [#l61101df] sageの式で変数として認識させるには変数をvar関数で宣言しなくてなりません。 変数の宣言は、 var('変数名') 複数の変数を宣言する場合には、スペースを空けて指定します。 宣言される変数を参照したりする場合には、 x, y = var('x y') のようにpython変数x, yに宣言したsage変数を代入します。 変数の値を指定して、関数の返す値を調べるには、 f1(x=1.5) のように式を表す変数にカッコを付けて、x=1.5と値を指定するとその時の関数の値が表示されます。 sage入力: #pre{{ x = var('x') f1 = (x - 1)*(x^2 - 1) f1(x=1.5) }} &color(blue){0.625000000000000}; 変数をx, yと複数存在する場合の例をです。f5に\(f5(x,y)=x2+y\)を定義し、f5(x=2, y=3)でx=2, y=3の時の値を出力します。 sage入力: #pre{{ var('x y w') f5=x^2 + y f5(x=2, y=3) }} &color(blue){7}; ** 規則の代入 [#u0590089] 規則の代入には、subs_exprを使います。 sage入力: #pre{{ f6=x^2 + 1; f6 }} &color(blue){x^2 + 1}; sage入力: #pre{{ f6.subs_expr(x^2 == w) }} &color(blue){w + 1}; ** 関数の定義 [#z365a917] sageで関数を定義する方法が2つあります。 - pythonの関数を定義する - lambda式を使って関数を定義する です。 *** python関数の定義 [#kdb65a09] 手続き的な関数を定義する場合には、python関数を使うと便利です。 pythonでは #pre{{ def 関数名(引数のリスト): 関数の本体 }} のように定義します。 以下にx3を返す関数cubeをpythonの関数を使って定義しました。 sage入力: #pre{{ def cube(x): return x^3 cube(2) }} &color(blue){8}; #pre{{ x, y = var('x y') cube(x+y) }} &color(blue){(x + y)^3}; *** lambda関数の定義 [#cb427343] つぎにlambda式を使って関数を定義する方法を示します。 x3のように式で表現できる場合には、lambda関数が便利です。 lambda 変数のリスト: 式 のように定義します。 先ほどと同じx3を返す関数の例を以下に示します。 sage入力: #pre{{ f = lambda x: x^3 f(2) }} &color(blue){8}; #pre{{ f(x+y) }} &color(blue){(x + y)^3}; ** 方程式の解法 [#a78e7648] 数式処理で助かる機能は、方程式の解法でしょう。 方程式の解法は、solve関数を使って以下のように行います。 solve(解放したい関数また関数のリスト, 解を得る変数) 例として、\(f(x)=ax+b\)の一次方程式を解いてみます。結果から\(x=−b/a\) が得られたので、\(a=2,b=1\)を代入して解の値を取得します。 sage入力: #pre{{ # 方程式の解法 var('x a b') f = a*x + b sln = solve(f, x); sln }} &color(blue){[x == -b/a]}; #pre{{ x = -b/a x(a=2, b=1) }} &color(blue){-1/2}; これだと、解から式をもう一度入力しなくてはなりません。そこで、a,bも変数として関数を定義し、solveでa=2,b=1をセットするようにします。こうすれば簡単に解の値が確認できます。 sage入力: #pre{{ var('x a b') g = lambda x, a, b : a*x + b solve(g(x, a, b), x) }} &color(blue){[x == -b/a]}; #pre{{ solve(g(x, a=2, b=1), x) }} &color(blue){[x == (-1/2)]}; 解の値がほしいときには、find_root関数で数値解析で値を計算します。 find_root(関数, 計算する変数の開始値、終了値) 先ほどの関数g(x):a=2,b=1の解をx=−1,1 の範囲で数値解析します。 sage入力: #pre{{ find_root(g(x, a=2, b=1), -1, 1) }} &color(blue){-0.5}; *** 多項式の解 [#u06bdee2] 一次方程式ではあまりに簡単なので、多項式の解を求めてみます。 $$ f(x)=x^3−2x+4 $$ のプロット図を以下に表示します。x軸と交差しているのはx=−2 です。-1と1付近に極値があります。 sage入力: #pre{{ var('x') f = x^3-2*x+4 plot(f, -2.5, 2.5) }} &ref(sage0.png); solveの結果、\(x=−2,x=1−i,x=1+i \)を得ました。 sage入力: #pre{{ solve(f, x) }} &color(blue){[x == (-I + 1), x == (I + 1), x == -2]}; sage入力: #pre{{ factor(f) }} &color(blue){(x + 2)*(x^2 - 2*x + 2)}; plotで可視化することでfind_rootでの計算範囲を大まかに把握することができます。 図から-3から3の間に解があることが見て取れたので、この範囲で数値解を求めます。 sage入力: #pre{{ find_root(f, -3, 3) }} &color(blue){-1.9999999999999949}; *** 連立方程式の解 [#q37df93e] 同様にsolveに関数に式のリストを与えることで、連立方程式の解を得ることができます。 $$ \begin{eqnarray} x^2 - 2y &=& 2 \\ x - y &=& 1 \end{eqnarray} $$ の解を例に以下に示します。 sage入力: #pre{{ var('x y') f = [y == 1/2*x^2 - 1, y == x - 1] solve(f, x, y) }} &color(blue){[ [x == 2, y == 1], [x == 0, y == -1] ]}; 上記関数のプロット結果は、以下のとおりです。 sage入力: #pre{{ p1 = plot(lambda x : 1/2*x^2 - 1, (x, -1, 3)) p2 = plot(lambda x : x - 1, (x, -1, 3)) (p1+p2).show() }} &ref(sage0-1.png); もう一つ、 $$ \begin{eqnarray} x^2 + y^2 &=& 1 \\ y &=& x - 1 \end{eqnarray} $$ の解を例に以下に示します。 sage入力: #pre{{ var('x y') f = [x^2 + y^2 == 1, y == x - 1] solve(f, x, y) }} &color(blue){[ [x == 1, y == 0], [x == 0, y == -1] ]}; この関数をimplicit_plotを使って関係式の形をそのまま使ってプロットしてみます。 円のなので、aspect_ratio=1としました。 sage入力: #pre{{ 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) }} &ref(sage0-2.png); ** 導関数 [#v361bc68] 高校で習った、導関数をsageを使って導いてみましょう。 関数\(f(x)=\frac{1}{2} x^3\)とします。 平均変化率は、\(g(x)=hf(x+h)−f(x)\)を以下のように定義します。 sage入力: #pre{{ # 導関数 def f(x): return (x^3)/2 var('x, h') g=(f(x + h) - f(x))/h; g }} &color(blue){1/2*((h + x)^3 - x^3)/h}; g(x)を展開し、整理すると sage入力: #pre{{ g1= g.rational_simplify(); g1 }} &color(blue){1/2*h^2 + 3/2*h*x + 3/2*x^2}; となります。 ここで、h→0 の極値を取ったものが求める導関数\(f′(x)\)です。 sage入力: #pre{{ limit(g1, h=0) }} &color(blue){3/2*x^2}; 結果は、f(x)をxで微分したものと同じになります。 sage入力: #pre{{ diff(f(x), x) }} &color(blue){3/2*x^2}; ** 微分 [#gc56b461] 高校で習う微分の公式は、 - \(\frac{d}{dx} c = 0\) - \(\frac{d}{dx} x^n = n x^{x-1}\) です。 sageの微分diff関数は、以下のように使います。 diff(関数, 微分する変数) 上記公式のsageで結果は、以下のようになります。 sage入力: #pre{{ # 微分 x, c, n = var('x c n') f = function('f', x) g = function('g', x) print diff(c, x) print diff(x^n, x) }} &color(blue){0}; &color(blue){n*x^(n - 1)}; sageで変数xの関数fを以下のように定義し微分すると、\(f′\)は以下のようにD[0](f)(x)となります。 #pre{{ x = var('x') f = function('f', x) }} これで変数xの関数fを定義したことになります。 sage入力: #pre{{ diff(f, x) }} &color(blue){D[0](f)(x)}; このD[0](f)を\(f′\)と読み替えると、 sage入力: #pre{{ 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) }} &color(blue){c*D[0](f)(x)}; &color(blue){D[0](f)(x) + D[0](g)(x)}; &color(blue){D[0](f)(x)*g(x) + f(x)*D[0](g)(x)}; &color(blue){(D[0](f)(x)*g(x) - D[0](g)(x)*f(x))/g(x)^2}; &color(blue){-D[0](g)(x)/g(x)^2}; &color(blue){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)\) *** いろんな関数の微分 [#l9d4ca27] 以下の関数を微分した結果を示します。 - \( (x^2+1)^3\) - \( sin(x)^3\) - \( sin^{-1}(x)\) - \( e^x \) - \( log(x)\) sage入力: #pre{{ # いろんな関数の微分 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) }} &color(blue){6*(x^2 + 1)^2*x}; &color(blue){3*sin(x)^2*cos(x)}; &color(blue){1/sqrt(-x^2 + 1)}; &color(blue){e^x}; &color(blue){1/x}; *** 微分の応用 [#ob786c66] 微分の応用として、接線の方程式とその法線を計算してみましょう。\( 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入力: #pre{{ f = exp(-(x^2)); f }} &color(blue){e^(-x^2)}; #pre{{ g = diff(f, x); g }} &color(blue){-2*x*e^(-x^2)}; #pre{{ m = g(x=1); m }} &color(blue){-2*e^(-1)}; #pre{{ 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) }} &ref(sage0-3.png); ** 積分 [#p3cac054] sageでは積分は、integrate関数で計算します。 #pre{{ integrate(被積分関数, 積分変数) また integrate(被積分関数, 積分変数, 積分範囲) }} とすると定積分を計算します。 \(\int xdx\)を計算した後、その微分を求めています。 sage入力: #pre{{ f = integrate(x); f }} &color(blue){1/2*x^2}; #pre{{ diff(f) }} &color(blue){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入力: #pre{{ # いろんな関数の積分 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) }} &color(blue){c*x}; &color(blue){log(x)}; &color(blue){e^x}; &color(blue){x*log(x) - x}; &color(blue){-cos(x)}; &color(blue){log(f(x))}; *** 定積分 [#oa2182e4] 次に定積分\( \int_0^{\pi/2}sin(x)dx \) の計算を示します。 sage入力: #pre{{ integrate(sin(x), x, 0, pi/2) }} &color(blue){1}; もう一つ $$ \int_0^3 x^2 + 2 sind(x) dx $$ の例を示し、 sage入力: #pre{{ integrate(x^2+2*sin(x), x, 0, 3) }} &color(blue){-2*cos(3) + 11}; 上記結果を数値とすると以下のようになります。 #pre{{ N(_) }} &color(blue){12.9799849932009}; *** 数値積分 [#l2f0a04e] 数値積分は、numerical_integral関数を使って計算します。 numerical_integral(被関分館数, 積分範囲) で戻り値は、積分結果と誤差を返します。 以下は、\(sin(x\))を0から\(\pi\)で数値積分した結果です。 sage入力: #pre{{ sigma, error = numerical_integral(sin(x), 0, pi); (sigma, error) }} &color(blue){(1.9999999999999998, 2.2204460492503128e-14)}; *** 積分の応用 [#hbbb8e1c] 積分の応用例として、サイクロイドの計算をしてみます。 サイクロイドは、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入力: #pre{{ t = var('t') x = 2*(t-sin(t)) y = 2*(1-cos(t)) parametric_plot([x, y], (t, 0, 2*pi)) }} &ref(sage0-4.png); #pre{{ cycloid = sqrt(diff(x,t)^2+diff(y,t)^2) }} 残念ながら、定積分は結果がでませんでした。 #pre{{ integrate(cycloid, t, 0, 2*pi) }} &color(blue){integrate(sqrt(4*(cos(t) - 1)^2 + 4*sin(t)^2), t, 0, 2*pi)}; 数値積分で、期待した結果となりました。 #pre{{ numerical_integral(cycloid, 0, 2*pi) }} &color(blue){(15.999999999999998, 1.7763568394002502e-13)}; ** コメント [#s902e4bb] #vote(おもしろかった[15],そうでもない[0],わかりずらい[1]) #vote(おもしろかった[16],そうでもない[0],わかりずらい[1]) 皆様のご意見、ご希望をお待ちしております。 #comment_kcaptcha