代数と微積分の基本#

Sageでは，初等的な代数と微積分に関連した多様な演算を実行することができる．例として，方程式の解を求める，微分や積分を計算する，ラプラス変換の実行などがあげられる． Sage Constructions には，さらに多様な具体例が盛られている．

オイラーによる連立微分方程式の解法#

次の例では，1階および2階微分方程式に対するオイラーの解法を具体的に解説する．手始めに1階微分方程式に対する解法の基本的アイデアを復習しておこう．初期値問題が

\[y'=f(x,y), \quad y(a)=c,\]

のような形式で与えられており， \(b>a\) を満足する \(x=b\) における解の近似値を求めたいものとする．

微分係数の定義から

\[y'(x) \approx \frac{y(x+h)-y(x)}{h},\]

ここで \(h>0\) は与えるべき小さな量である．この近似式と先の微分方程式を組み合わせると \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\) が得られる．これを \(y(x+h)\) について解くと:

\[y(x+h) \approx y(x) + h\cdot f(x,y(x)).\]

(他にうまい呼び方も思いつかないので) \(h \cdot f(x,y(x))\) を "補正項" と呼び， \(y(x)\) を \(y\) の "更新前項(old)", \(y(x+h)\) を \(y\) の "更新後項(new)"と呼ぶことにすると，上の近似式を

\[y_{new} \approx y_{old} + h\cdot f(x,y_{old}).\]

と表わすことができる．

ここで \(a\) から \(b\) までの区間を \(n\) ステップに分割すると \(h=\frac{b-a}{n}\) と書けるから，ここまでの作業から得られた情報を整理して以下の表のようにまとめることができる．

\(x\)	\(y\)	\(h\cdot f(x,y)\)
\(a\)	\(c\)	\(h\cdot f(a,c)\)
\(a+h\)	\(c+h\cdot f(a,c)\)	...
\(a+2h\)	...
...
\(b=a+nh\)	???	...

我々の目標は，この表の空欄を上から一行づつ全て埋めていき，最終的に \(y(b)\) のオイラー法による近似である???に到達することである．

連立微分方程式に対する解法もアイデアは似ている．

例題: \(z''+tz'+z=0\), \(z(0)=1\), \(z'(0)=0\) を満足する \(t=1\) における \(z(t)\) を，4ステップのオイラー法を使って数値的に近似してみよう．

ここでは問題の2階常微分方程式を( \(x=z\), \(y=z'\) として)二つの1階微分方程式に分解してからオイラー法を適用することになる。

Sage

sage: t,x,y = PolynomialRing(RealField(10),3,"txy").gens()
sage: f = y; g = -x - y * t
sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
      t                x            h*f(t,x,y)                y       h*g(t,x,y)
      0                1                  0.00                0           -0.25
    1/4              1.0                -0.062            -0.25           -0.23
    1/2             0.94                 -0.12            -0.48           -0.17
    3/4             0.82                 -0.16            -0.66          -0.081
      1             0.65                 -0.18            -0.74           0.022

Python

>>> from sage.all import *
>>> t,x,y = PolynomialRing(RealField(Integer(10)),Integer(3),"txy").gens()
>>> f = y; g = -x - y * t
>>> eulers_method_2x2(f,g, Integer(0), Integer(1), Integer(0), Integer(1)/Integer(4), Integer(1))
      t                x            h*f(t,x,y)                y       h*g(t,x,y)
      0                1                  0.00                0           -0.25
    1/4              1.0                -0.062            -0.25           -0.23
    1/2             0.94                 -0.12            -0.48           -0.17
    3/4             0.82                 -0.16            -0.66          -0.081
      1             0.65                 -0.18            -0.74           0.022

したがって， \(z(1)\approx 0.65\) が判る．

点 \((x,y)\) をプロットすれば、その曲線としての概形を見ることができる．それには関数 eulers_method_2x2_plot を使うが，その前に三つの成分(\(t\), \(x\), \(y\))からなる引数を持つ関数 \(f\) と \(g\) を定義しておかなければならない．

Sage

sage: f = lambda z: z[2]        # f(t,x,y) = y
sage: g = lambda z: -sin(z[1])  # g(t,x,y) = -sin(x)
sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)

Python

>>> from sage.all import *
>>> f = lambda z: z[Integer(2)]        # f(t,x,y) = y
>>> g = lambda z: -sin(z[Integer(1)])  # g(t,x,y) = -sin(x)
>>> P = eulers_method_2x2_plot(f,g, RealNumber('0.0'), RealNumber('0.75'), RealNumber('0.0'), RealNumber('0.1'), RealNumber('1.0'))

この時点で， P は2系列のプロットを保持していることになる． \(x\) と \(t\) のプロットである P[0] ，および \(y\) と \(t\) のプロットである P[1] である．これら二つをプロットするには、次のようにする:

Sage

sage: show(P[0] + P[1])

Python

>>> from sage.all import *
>>> show(P[Integer(0)] + P[Integer(1)])

(プロットの詳細についてはプロットする節を参照．)

特殊関数#

数種類の直交多項式と特殊関数が，PARI [GAP] およびMaxima [Max] を援用して実装されている．詳細についてはSageレファレンスマニュアルの“Orthogonal polynomials"(直交多項式)と“Special functions"(特殊関数)を参照してほしい．

Sage

sage: x = polygen(QQ, 'x')
sage: chebyshev_U(2,x)
4*x^2 - 1
sage: bessel_I(1,1).n(250)
0.56515910399248502720769602760986330732889962162109200948029448947925564096
sage: bessel_I(1,1).n()
0.565159103992485
sage: bessel_I(2,1.1).n()
0.167089499251049

Python

>>> from sage.all import *
>>> x = polygen(QQ, 'x')
>>> chebyshev_U(Integer(2),x)
4*x^2 - 1
>>> bessel_I(Integer(1),Integer(1)).n(Integer(250))
0.56515910399248502720769602760986330732889962162109200948029448947925564096
>>> bessel_I(Integer(1),Integer(1)).n()
0.565159103992485
>>> bessel_I(Integer(2),RealNumber('1.1')).n()
0.167089499251049

ここで注意したいのは，Sageではこれらの関数群が専ら数値計算に便利なようにラップ(wrap)されている点だ．記号処理をする場合には，以下の例のようにMaximaインターフェイスをじかに呼び出してほしい．

Sage

sage: maxima.eval("f:bessel_y(v, w)")
'bessel_y(v,w)'
sage: maxima.eval("diff(f,w)")
'(bessel_y(v-1,w)-bessel_y(v+1,w))/2'

Python

>>> from sage.all import *
>>> maxima.eval("f:bessel_y(v, w)")
'bessel_y(v,w)'
>>> maxima.eval("diff(f,w)")
'(bessel_y(v-1,w)-bessel_y(v+1,w))/2'

代数と微積分の基本#

方程式を解く#

方程式を解析的に解く#

方程式を数値的に解く#

微分，積分，その他#

微分方程式を解く#

オイラーによる連立微分方程式の解法#

特殊関数#