計算機基礎A - 14・15: Maxima の使い方
Table of Contents
1 Maxima とは
Maxima （マキシマ）はフリーで使える数式処理システムです． 微分積分や線形代数の計算だけでなく，式変形やテイラー展開なども簡単に行うことができます． 数学の勉強やレポート作成などに威力を発揮するでしょう． Maxima の情報は以下の Web ページから得ることができます．
他にもインターネット上に多くの情報がありますので，検索してみるとよいでしょう．
Ubuntu や Debian を利用している場合は apt-get で簡単にインストールできます．
ターミナルから次のようにコマンドを入力して下さい．
sudo apt-get install wxmaxima maxima-emacs
2 起動と終了
2.1 起動方法
Maxima をターミナルから使ってみましょう．
ターミナルから maxima というコマンドを実行すると Maxima が起動します．
$ maxima Maxima 5.20.1 http://maxima.sourceforge.net using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (a.k.a. GCL) Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1)
(%i1) というのは Maxima のプロンプトです．
Maxima では入力と出力に通し番号が付けられます．
N 番目の入力は %iN ，出力は %oN となります．
ユーザはこのラベルによって以前の入力と出力を参照（再利用）することができます．
試しに 2^100 （2 の 100 乗）を計算させてみましょう．
Maxima のプロンプトに対して 2^100; と入力し，リターンキー（エンターキー）を押します．
(%i1) 2^100; (%o1) 1267650600228229401496703205376
最後の ; が入力の区切りになります．
これを忘れると入力が継続されていると判断され，結果が出力されませんので，注意しましょう．
では，さきほどの結果が正しいか検算してみましょう．
1番目の入力 %i1 に対する出力結果は %o1 で参照することができます．
それを 1/100 乗して 2 に戻るかチェックすればよいですね．
(%i2) %o1^(1/100); (%o2) 2
Maxima を使っていると，直前の計算結果を参照したいことがよくあります．
これは番号なしの % で参照することができます．
(%i3) % + x; (%o3) x + 2
% が 2 に置き換えられたことに注意してください．
また x という記号は記号のまま計算されています．
このように，Maxima では記号を使った計算（式の展開や因数分解，不定積分，テイラー展開など）をすることもできます。
以下は多項式の計算例です．
(%i4) % + (2*x+3); (%o4) 3 x + 5
2.2 終了方法
Maxima を終了するには quit() コマンドを入力します．
Maxima が終了すると，シェルのプロンプトが戻ってきます．
(%i5) quit();
$
2.3 Maxima を Emacs から使う
Maxima は Emacs の maxima-mode から使うことができます．
Emacs の編集機能が利用できるので，ターミナルから使うよりずっと便利です．
Emacs から M-x maxima RET と入力してみましょう．
RET はリータンキーまたはエンターキーを押すことを示します．
2.4 GUI な Maxima
最近は Mathematica や Maple といった市販の数式処理システムに似た
GUI を備えた Maxima を利用することもできるようです．
ターミナルから wxmaxima というコマンドを実行すると GUI 付きの Maxima が起動します．
数式がきれいに表示されますので，試してみるとよいでしょう．
3 使ってみよう
3.1 数式の記号
Maxima には四則演算に加えいくつかの記号が用意されています．
+ 加算 - 減算 * 乗算 / 算除 ^ ベキ乗 ! 階乗 !! 二重階乗
加減乗除 +, -, *, / の優先順序はふつうの数学と同じです．
優先順序を変更したい場合は丸括弧 ( ... ) でグルーピングします．
(%i1) 1 + 2 * 3; (%o1) 7 (%i2) (1 + 2) * 3; (%o2) 9
x^2 は x の二乗を表します．
2 の 10 乗を計算するには次のように入力します．
(%i3) 2^10; (%o3) 1024
n!! は二重階乗（1つおきの階乗）を計算します．
(%i4) 10!!; /* 2*4*6*8*10 */ (%o4) 3840 (%i5) 11!!; /* 1*3*5*7*9*11 */ (%o5) 10395
3.2 代入
変数に値を代入するにはセミコロン : を使います．
変数名: 値
変数 a に 2/3 を代入するには次のように入力します．
(%i1) a: 2/3; 2 (%o1) - 3
変数に値を代入しておくと，あとでその変数を利用することができます．
(%i2) a - 1/6; 1 (%o2) - 2
変数を単なる記号に戻したい場合など，変数から値を消去するには remvalue() を使います．
(%i3) remvalue(a);
(%o3) [a]
すべての変数を消去するには kill(all) とします．
4 数の計算
4.1 代数計算
Maxima は加減乗除やベキ乗の計算をすることができます．
(%i1) 123 + 456; (%o1) 579 (%i2) 123 * 456; (%o2) 56088
整数を扱っているとき Maxima は厳密な値を返します． 電卓では扱えないような大きな数でも大丈夫です．
(%i3) 2^456; (%o3) 186070713419675363980626894819329160794532188335953423432061490990243657\ 757029868371504908982723472783555205531204141550984858016925351936 (%i4) 100!; (%o4) 933262154439441526816992388562667004907159682643816214685929638952175999\ 932299156089414639761565182862536979208272237582511852109168640000000000000000\ 00000000
4.2 有理数
2つの整数の商として表される数を「有理数」といいます． 有理数のみを含む計算では，結果を有理数または整数として返します．
(%i1) 1/2 + 1/3; 5 (%o1) - 6 (%i2) % + 1/6; (%o2) 1
1/5 という有理数は数学的には 0.2 に等しいですが，Maxima はこれらを区別して扱います． すなわち，整数と有理数は何桁でも厳密に扱いますが，小数点のついた数（浮動小数点数）は「近似値」として扱います． この区別は重要ですから，しっかり覚えておきましょう．
(%i3) 1/5 + 1; 6 (%o3) - 5 (%i4) 0.2 + 1; (%o4) 1.2
4.3 無理数
Maxima は結果をなるべく厳密に保とうとするので，無理数は無理数のまま結果を返します．
例えば，3 の平方根 √(3) を求めると，Maxima は 1.732… ではなく √(3) を表す式
sqrt(3) を返します．
(%i1) sqrt(3); (%o1) sqrt(3)
一方，近似値である浮動小数点数 3.0 の平方根 sqrt(3.0) を求めると，結果も浮動小数点数で返されます．
(%i2) sqrt(3.0); (%o2) 1.732050807568877
4.4 近似値としての小数
Maxima では浮動小数点数は近似値とみなされます． 式の中に浮動小数点数が含まれていれば，Maxima はその式全体を近似値として扱い，結果を浮動小数点数で表します． 次の例を見てください．
(%i1) 1/2 + 3.4/5; (%o1) 1.18
この例では， 1/2 は厳密な数ですが 3.4 は近似値です．
そのため，実際の計算は次のように行われます．
3.4/5が0.68と計算される1/2 + 0.68が0.5 + 0.68に変換される （有理数1/2が浮動小数点数0.5に変換される）1.18が結果として返される
浮動小数点数（近似値）を整数や有理数（厳密な数）に変換するには次のような関数を使います．
ceiling(x) x 以上の最小の整数を返す floor(x) x 以下の最大の整数を返す round(x) x にもっとも近い整数を返す
4.5 複素数
Maxima は複素数を扱うことができます．
虚数単位 i は %i で表します．
(%i1) (1 + 2*%i) + (3 + 4*%i); (%o1) 6 %i + 4 (%i2) (1 + 2*%i) / (3 + 4*%i); 2 %i + 1 (%o2) -------- 4 %i + 3
分母を実数化するには algebraic を true にして ratsimp() を実行します．
最初は意味がわからないかもしれませんが，おまじないだと思ってください．
(%i3) ratsimp(%), algebraic: ture; 2 %i + 11 (%o3) --------- 25
複素数 z は2つの実数 x, y を用いて， z = x + i y と表すことができます． x のことを実部 (real part)，y のことを虚部 (imaginary part) といいます． 複素数 z の実部と虚部を求めるには次の関数を利用します．
realpart(z) 実部を返す imagpart(z) 虚部を返す
4.6 近似値を求める
Maxima は整数と有理数を厳密に扱いますが，関数 float(), bfloat()
を使うと近似値を求めることができます．
(%i1) float(sqrt(3)); (%o1) 1.732050807568877
float() は16桁の近似値を返します．
それ以上の桁数が必要な場合は bfloat() を使います．
bfloat() の返す桁数は変数 fpprec で制御することができます．
3 の平方根を50桁ほど求めてみましょう．
(%i2) bfloat(sqrt(3)), fpprec: 50; (%o2) 1.7320508075688772935274463415058723669428052538104b0
4.7 定数
Maxima には次のような定数があらかじめ定義されています．
%pi 円周率 %e 自然対数の底 %i 虚数単位 inf 正の無限大
これらの記号は数学的に厳密な値を表します． 特に指定しない限り，Maxima はこれらの値の近似値を計算しません．
(%i1) %pi * 2; (%o1) 2 %pi (%i2) bfloat(%pi), fpprec: 100; (%o2) 3.1415926535897932384626433832795028841971693993751058209749445923078164\ 06286208998628034825342117068b0
4.8 基本的な関数
Maxima には次のような関数が用意されています．
sqrt(x) 平方根 exp(x), log(x) 指数関数・対数関数 sin(x), cos(x), tan(x) 三角関数 asin(x), acos(x), atan(x) 三角関数の逆関数 sinh(x), cosh(x), tanh(x) 双曲線関数 abs(x) 絶対値
sin(%pi/2) の値を求めてみましょう．
(%i1) sin(%pi/2); (%o1) 1
4.9 問題
- 365日は何秒か？
bfloat()を用いて \(\pi\) の値を 300 桁まで求めよ．e^π(\(e^{\pi}\)) とπ^e(\(\pi^{e}\)) はどちらが大きいか？
5 式を使った計算
5.1 数式の計算
1 から 10 までの和を求めてみましょう．
(%i1) 1+2+3+4+5+6+7+8+9+10; (%o1) 55
これは公式 1+2+3+...+n = n*(n+1)/2 を使って計算することができます．
(%i2) n*(n+1)/2, n=10; (%o2) 55
n*(n+1)/2, n=10 は ev(n*(n+1)/2, n=10) の省略形で，式に一時的に値を代入することができます．
ev は evaluate （評価する）のことでしょう．
5.2 数列の和（総和）
数列の和を sum(), nusum() で計算できます．
sum(一般項, 添字変数, 下限, 上限) nusum(一般項, 添字変数, 下限, 上限)
1 + 2 + ... + 10 を今度は sum() を使って計算してみましょう．
添字変数を i とすると，一般項は i ，添字の下限は 1 ，上限 10 となります．
(%i12) sum(i, i, 1, 10); (%o12) 55
同様に， 1^2 + 2^2 + ... + 10^2 を計算してみましょう．
(%i13) sum(i^2, i, 1, 10); (%o13) 385
nusum() を使うと第 n 項（n は記号）までの和を計算することもできます．
(%i14) nusum(i, i, 1, n); n (n + 1) (%o14) --------- 2 (%i15) nusum(i^2, i, 1, n); n (n + 1) (2 n + 1) (%o15) ------------------- 6 (%i16) nusum(i^3, i, 1, n); 2 2 n (n + 1) (%o16) ----------- 4
5.3 数列の積（総積）
Maxima には総積 Π （大文字のパイ）を計算する関数も用意されています．
product(式, 添字, 下限, 上限)
product() を使って sin(1*π/4) * sin(2*π/4) * sin(3*π/4) を計算してみましょう。
(%i5) product(sin(k*%pi/4), k, 1, 3); 1 (%o5) - 2
5.4 積分
Maxima は積分を計算することができます．
関数 integrate() は f(x) を x について min から max まで積分します．
integrate(f(x), x, min, max)
sin(x) を区間 [0, π] で積分してみましょう．
(%i1) integrate(sin(x), x, 0, %pi); (%o1) 2
積分についてはまたあとで説明します．
5.5 方程式の近似解
Maxima には方程式の根（方程式の解）を数値的に求める関数が用意されています． 答えは浮動小数点数（近似値）になります． （あとで厳密解を求める方法も紹介します．）
allroots(f(x)) 実数係数の方程式 f(x)=0 の実数解と複素解を求める find_root(f(x), x, min, max) f(x)=0 の解を区間 [a, b] の間でみつける
まず， allroots() を利用して二次方程式 x^2 - x - 2 = 0 の解を求めてみましょう．
(%i1) allroots(x^2-x-2); (%o1) [x = - 1.0, x = 2.0]
解くべき方程式は f(x) = a という形で与えることもできます．
(%i2) allroots(x^2 = 4); (%o2) [x = 2.0, x = - 2.0]
allroots() はすべての解を探しますが，
find_root() を使うとある区間に含まれる解だけを探すことができます．
x^2 - x - 2 = 0 の解を 0 から 5 の範囲で探してみましょう．
(%i3) find_root(x^2-x-2, x, 0, 5); (%o3) 2.0
5.6 問題
sum()を用いて1 + 3 + 5 + ... + 99を計算せよ。sum()を用いて1^3 + 2^3 + ... + 10^3を計算せよ．sum()を用いて1/1 + 1/2! + 1/3! + ... + 1/10!を計算せよ．product()を用いて 10! （つまり 1*2*…*10）を求めよ．exp(1*π/6) * exp(2*π/6) * exp(3*π/6) * exp(4*π/6) * exp(5*π/6)を計算せよ．- 関数
1/(1+x^2)を区間[0, 1]で積分せよ． x^3 - 4*x^ = 0の解を求めよ．- 方程式
sin(x) = cos(x)の解を求めよ．（ヒント：区間[0, %pi]で探せ） - 方程式
x * tan(x) = 1の解を求めよ．（ヒント：区間[0, %pi/2]で探せ）
6 記号計算
Maxima では x や y など，記号を含んだ式をそのまま扱うことができます． 式などを操作する関数には次のようなものがあります．
expand() 式の展開 factor() 因数分解 partfrac() 部分分数展開 ratsimp() 式の簡単化
6.1 展開と因数分解
(x+y)^4 を展開してみましょう．
(%i1) expand((x+y)^4); 4 3 2 2 3 4 (%o1) y + 4 x y + 6 x y + 4 x y + x
結果を因数分解してみましょう． ちゃんと元に戻るでしょうか？
(%i2) factor(%); 4 (%o2) (y + x)
因数が違っていても因数分解してくれます．
x^2 - 3*x + 2 で確かめてみましょう．
(%i3) factor(x^2 - 3*x + 2); (%o3) (x - 2) (x - 1)
6.2 部分分数展開・式の簡単化
x/(x+1) を部分分数展開してみましょう。
(%i4) partfrac(x/(x+1), x); 1 (%o4) 1 - ----- x + 1
ratsimp() でこの式を簡単化します。
(%i5) ratsimp(%); x (%o5) ----- x + 1
6.3 問題
x^2 + 5*x + 6を因数分解せよ．(x-y)^2 - 6*(x-y) + 8を因数分解せよ．a = 1/(sqrt(7) - sqrt(5)),b = 1/(sqrt(7) + sqrt(5))のとき， a と b の積a * bの値を求めよ． ヒント：ratsimp(a*b)で式を簡単化すればよい．
7 方程式を厳密に解く
allroots() と find_root() で方程式の近似解を求めましたが，
solve() を利用すると方程式を厳密に解くこともできます．
solve(方程式, 変数) solve([方程式1, 方程式2, ...], [変数1, 変数2, ...])
7.1 方程式の解
solve() を利用して，方程式 x^2 - 1 = 0 を解いてみましょう．
(%i1) solve(x^2 - 1 = 0, x); (%o1) [x = - 1, x = 1]
7.2 連立方程式
連立方程式を解くこともできます．
連立方程式の場合は，連立させる方程式を角括弧 [ ... ] の中にカンマで区切って並べます．
y = x^2, y = x+2 を解いてみましょう．
(%i2) solve([y=x^2, y=x+2], [x, y]); (%o2) [[x = - 1, y = 1], [x = 2, y = 4]]
7.3 解けない場合
ただし， solve() は万能ではありません．
例えば， x * exp(x) = 1 という方程式を考えてみましょう．
解が区間 [0, 1] に含まれることは明らかですが，それを solve() で求めることはできません．
このような場合には find_root() で答え（近似解）を見つけましょう．
(%i4) solve(x*exp(x) = 1, x); - x (%o4) [x = %e ] (%i5) find_root(x*exp(x) = 1, x, 0, 1); (%o5) 0.56714329040978
7.4 問題
- 方程式
x^2 + 1 = 0を解け． - 円
x^2 + y^2 = 10と直線y = 2*x-5の共有点の座標を求めよ．
8 関数の定義
8.1 関数定義の例
関数の定義には := を使います。
f(x) = x^2 + 2*x + 1 という二次関数を定義するには次のようにします．
(%i1) f(x) := x^2 + 2*x + 1; 2 (%o1) f(x) := x + 2 x + 1
8.2 関数の利用
これで f(x) は x^2 + 2*x + 1 とまったく同じ働きをするようになりました．
組み込み関数 factor() で因数分解してみましょう．
(%i2) factor(f(x)); 2 (%o2) (x + 1)
関数の括弧の中身（例えば f(x) の x） のことを「引数（ひきすう）」といいます．
関数を定義しておくと，引数に値を代入して，関数値を計算することができます．
(%i3) f(1/2); 9 (%o3) - 4
8.3 if 文による場合分け
場合分けが必要な場合は if 文を利用します． if 文には次のような書き方があります．
if 条件式 then 条件式が成り立つ場合の処理 if 条件式 then 条件式が成り立つ場合の処理 else 成り立たない場合の処理 if 条件式1 then 処理1 elseif 条件式2 then 処理2 ... elseif 条件式n then 処理n else 処理0
x < 0 で 0, x >= 0 で 1 となる単位ステップ関数 s(x) を定義してみましょう．
(%o1) s(x) := if x < 0 then 0 else 1; (%o1) s(x) := if x < 0 then 0 else 1
x < 0 で 0, x = 0 で 1/2, x > 0 で 1 となる関数 h(x) を定義するには次のようにします．
これはヘヴィサイド関数と呼ばれます．
(%i2) h(x) := if x < 0 then 0 elseif x = 0 then 1/2 else 1; 1 (%o2) h(x) := if x < 0 then 0 elseif x = 0 then - else 1 2
9 グラフの作成
9.1 2次元グラフ
関数 y = f(x) のグラフを作図するには plot2d() を利用します．
plot2d(f(x), [x, x の下限, x の上限]) plot2d([f1(x), f2(x), ..., fn(x)], [x, x の下限, x の上限]) (複数のグラフを重ねて表示する場合)
関数 y = sin(x)/x のグラフを [-10, 10] の範囲で表示してみましょう．
(%i1) plot2d(sin(x)/x, [x, -10, 10]); plot2d: expression evaluates to non-numeric value somewhere in plotting range. (%o1) /Users/takagi/maxout.gnuplot_pipes
こんどは関数 y = x と y = x * sin(x) のグラフを重ねて表示してみましょう．
(%i2) plot2d([x, x*sin(x)], [x, -10, 10]); (%o2) /Users/takagi/maxout.gnuplot_pipes
9.2 媒介変数表示された2次元グラフ
x, y 座標の値をあるパラメータを使って表すことを「媒介変数表示」といいます．
例えば， x = cos(t), y = sin(t) は半径 1 の円のパラメータ t による媒介変数表示です．
plot2d() は媒介変数表示されたグラフを表示することもできます．
plot2d([parametric, x(t), y(t)], [t, t の下限, t の上限], [nticks, グラフの分割数])
実際に，媒介変数表示された円のグラフを表示してみましょう．
(%i1) plot2d([parametric, cos(t), sin(t)], [t, 0, 2*%pi], [nticks, 360]); (%o1) /Users/takagi/maxout.gnuplot_pipes
9.3 陰関数によって表された2次元グラフ
x と y の関係が y について解かれておらず，
f(x, y) = 定数 の形で表されたものを陰関数といいます．
例えば，半径 1 の円の方程式 x^2 + y^2 = 1 は陰関数の一種です．
implicit_plot パッケージを利用すると，陰関数のグラフをプロットすることができます．
implicit_plot(f(x,y) = 定数, [x, x の下限, x の上限], [y, y の下限, y の上限])
円 x^2 + y^2 = 1 を図示してみましょう．
implicit_plot() を利用するには implicit_plot パッケージを読み込む必要があります．
パッケージの読み込みには load() コマンドを使います．
(%i3) load(implicit_plot)$ (%i4) implicit_plot(x^2 + y^2 = 1, [x, -1.2, 1.2], [y, -1.2, 1.2])$
この例のように，入力の末尾を ; ではなく $ にすると Maxima からの出力が抑制されます．
戻り値が不要なときは $ を使うとよいでしょう．
9.4 3次元グラフ
z = f(x, y) は三次元空間内の曲面を表す方程式です．
plot3d() を利用すると，簡単に曲面のグラフを図示することができます．
plot3d(f(x,y), [x, x の下限, x の上限], [y, y の下限, y の上限], [grid, x の分割数, y の分割数])
分割数を増やすことで滑らかな描画を得ることができます．
例として f(x,y) = sin(sqrt(x^2+y^2)) / sqrt(x^2+y^2) で表される曲面を図示してみましょう．
x と y の範囲はどちらも [-3π, 3π] とします．
以下では r(x,y) = sqrt(x^2 + y^2) という関数を定義して利用しています．
(%i1) r(x,y) := sqrt(x^2 + y^2); 2 2 (%o1) r(x, y) := sqrt(x + y ) (%i2) plot3d(sin(r(x,y))/r(x,y), [x, -3*%pi, 3*%pi], [y, -3*%pi, 3*%pi], [grid, 100, 100]); (%o2) /Users/takagi/maxout.gnuplot_pipes
9.5 練習
x(t) = sin(4*t)*cos(t),y(t) = sin(4*t)*sin(t)のように媒介変数表示されたグラフを[0, 2π]の範囲で表示せよ．x^2 - y^2 = 1のグラフを図示せよ． x と y の範囲はどちらも[-5, 5]とすること．
10 代数式を操作する
代数式を操作する関数として次のようなものが用意されています．
coeff(式, x, n) 式中の x^n の係数を返す num(有理式) 分子を返す denom(有理式) 分母を返す ratsubst(a, b, c) c 中の b に a を代入し，結果の式を返す gcd(式1, 式2, ...) 最大公約数を返す（式は多項式でもよい）
10.1 最大公約数
2つの整式 x^2+2*x+1, x^2+3*x+2 の最大公約数 (GCM) を求めてみましょう．
(%i1) gcd(x^2+2*x+1, x^2+3*x+2); (%o1) x + 1
上で得られた x + 1 の x に 9 を代入してみましょう．
(%i2) ratsubst(9, x, %); (%o2) 10
10.2 問題
gcd()を利用して，最小公倍数を返す関数lcm()を定義せよ．
11 三角関数を含む式の簡単化
三角関数を含む式の簡単化には次のような関数を利用します．
trigsimp() 三角関数・指数関数を含む式を簡単化する trigreduce() 三角関数同士の積を減らして簡単化する trigexpand() 加法定理や倍角公式を使って式を展開する
11.1 三角関数の簡単化の例
cos(x)^2 - sin(x)^2 を簡単化してみましょう．
（何を基準に「簡単」というか，難しいところですが）
(%i1) trigsimp(cos(x)^2 - sin(x)^2); 2 (%o1) 2 cos (x) - 1 (%i2) trigreduce(cos(x)^2 - sin(x)^2); cos(2 x) + 1 cos(2 x) 1 (%o2) ------------ + -------- - - 2 2 2 (%i3) trigsimp(%); (%o3) cos(2 x)
今度は cos(3*x) を展開してみましょう．
(%i4) trigexpand(cos(3*x)); 3 2 (%o4) cos (x) - 3 cos(x) sin (x) (%i5) trigsimp(%); 3 (%o5) 4 cos (x) - 3 cos(x)
11.2 問題
- 三角関数の加法定理
trigexpand(sin(x+y))などを確認せよ． - 双曲線関数の加法定理
trigexpand(sinh(x+y))などを確認せよ．
12 微分
関数や式を微分するには diff() を使います．
diff(f(x), x) f(x) を x で微分 diff(f(x), x, n) f(x) を x で n 階微分
12.1 微分の例
x^2 を微分してみましょう．
(%i1) diff(x^2, x); (%o1) 2 x
cos(x) や exp(x) などの関数を微分することもできます。
例として， x*cos(x) を微分してみましょう。
(%i2) diff(x*cos(x), x); (%o2) cos(x) - x sin(x)
なんと x^x （x の x 乗）を微分することもできます．
(%i3) diff(x^x, x); x (%o3) x (log(x) + 1)
これは正しい答えになっているのですが，わかりますか？
（ヒント： x^x = exp(x*log(x)) のように変形してから微分します）
12.2 問題
y = (x+2)*(3*x-1)を2階微分せよ．y = x * sqrt(x)を微分せよ．y = e^(x^2)を微分せよ． ヒント：%e^(x^2)またはexp(x^2)を微分する．
13 積分
integrate() を利用すると不定積分と定積分が計算できます．
不定積分の場合，積分定数は表示されないので注意してください．
integrate(f(x), x) 不定積分 integrate(f(x), x, a, b) 区間 [a, b] の定積分
また integrate() では積分がうまく計算できないこともあります．
そのような場合は，後述の romberg() を用いた数値積分を検討するとよいでしょう．
13.1 不定積分
f(x) = x の不定積分を求めてみましょう．
(%i5) integrate(x, x); 2 x (%o5) -- 2
正しい結果は x^2/2 + C （C は積分定数）です．
というわけで，定数を除けば正解ですね．
こんどは f(x) = x^2*cos(x) + 2*x*sin(x) の不定積分を求めてみましょう．
(%i6) integrate(x^2*cos(x) + 2*x*sin(x), x); 2 (%o6) 2 (sin(x) - x cos(x)) + (x - 2) sin(x) + 2 x cos(x)
この結果はもう少し簡単になりそうです．
このような場合は ratsimp() で式を簡単化できます。
(%i7) ratsimp(%o18); 2 (%o7) x sin(x)
次の積分は覚えておくと便利です．
(%i8) integrate(log(x),x); (%o8) x log(x) - x
13.2 定積分
関数 f(x) = sin(x) を 0 から π まで積分してみましょう．
円周率 π は %pi という定数で用意されています．
(%i8) integrate(sin(x), x, 0, %pi); (%o8) 2
inf という記号を使えば，積分の上限または下限を無限大にとることもできます．
次の積分はガウス積分と呼ばれ，物理では非常に重要です．
(%i9) integrate(exp(-x^2), x, -inf, inf); (%o9) sqrt(%pi)
13.3 問題
x+1の不定積分を求めよ．x+1を[0, 1]の範囲で積分せよ．
14 数値積分
解析的に積分できない場合でも romberg() によって数値積分を実行することができます．
romberg(f(x), x, 下限, 上限) 区間 [a, b] の定積分
14.1 数値積分の例
f(x) = sin(sin(x)) を 0 から 1 まで積分してみましょう．
次に示すように integrate() は積分に失敗しますが romberg() なら積分値が求められます．
(%i1) integrate(sin(sin(x)), x, 0, 1); 1 / [ (%o1) I sin(sin(x)) dx ] / 0 (%i2) romberg(sin(sin(x)), x, 0, 1); (%o2) 0.43060592364256
14.2 問題
cos(x + cos(x))を[0, 1]の範囲で積分せよ。exp(-x^2)を[-10, 10]の範囲で積分し，結果の2乗を求めよ．
15 極限
limit() を用いて極限値を計算することができます．
limit(式, 変数, 値) limit(式, 変数, 値, 方向)
方向 には極限を取る向きを plus または minus で指定します．
上からの極限が plus ，下からの極限が minus です．
15.1 極限の例
x を 0 に近づけたときの sin(x)/x 極限値を求めてみましょう．
(%i1) limit(sin(x)/x, x, 0); (%o1) 1
無限大での極限を計算することもできます．
(%i2) limit(log(x)/x, x, inf); (%o2) 0
15.2 問題
- x を
正の無限大に近づけたときの(1+1/x)^xの極限値を求めよ． - x を
+0に近づけたときのx*log(x)の極限値を求めよ．+0は上から（つまり，正の方向から） 0 に近づくことを表します．
16 テイラー展開
taylor() を使うと関数のテイラー展開を求めることができます．
taylor(f(x), x, a, n) f(x) を x=a の周りで n 次の項までテイラー展開する
16.1 テイラー展開の例
sin(x) を x=0 の周りで5次の項まで展開してみましょう．
(%i1) taylor(sin(x), x, 0, 5); 3 5 x x (%o1)/T/ x - -- + --- + . . . 6 120
coeff() を使って3次の項の係数を求めてみましょう．
(%i2) coeff(%, x, 3); 1 (%o2)/R/ - - 6
16.2 問題
exp(x)をx=0の周りで5次の項までテイラー展開せよ．cos(x)をx=0の周りで5次の項までテイラー展開せよ．
17 ベクトル
ベクトルを定義するには角括弧 [ ... ] を使います．
[要素1, 要素2, ..., 要素n]
17.1 ベクトルの例
次の例は，変数 u に (1, 2, 3) という三次元ベクトルを代入しています．
(%i1) u: [1, 2, 3]; (%o1) [1, 2, 3]
ベクトル u の絶対値を求めてみましょう．
ベクトル u, v の内積は u . v で計算できるので，
自分自身との内積 u . u はベクトル u の絶対値の二乗になることを利用しましょう．
(%i2) u . u; (%o2) 14 (%i3) sqrt(%); (%o3) sqrt(14)
したがって、ベクトルの絶対値を返す関数 norm() は次のように定義できます．
(%i4) norm(x) := sqrt(x . x); (%o4) norm(x) := sqrt(x . x) (%i5) norm(u); (%o5) sqrt(14)
17.2 問題
- ベクトル
a = [1, 1, -1],b = [2, -1, 1]の内積を求めよ． - ベクトル
a = [1, 1, 1],b = [1, -1, 1]のなす角を θ とする．cosθの値を求めよ．
18 行列
18.1 行列の定義
行列は matrix() で定義します．
matrix([1行目の成分リスト], [2行目の成分リスト], ..., [n行目の成分リスト])
例として，次のような2行2列の正方行列 A を考えます．
A = ( 1 -2 )
(-3 4 )
この行列を定義するには次のようにします． 要素を並べる順番に注意しましょう．
(%i1) A: matrix([1, -2], [-3, 4]); [ 1 - 2 ] (%o1) [ ] [ - 3 4 ]
18.2 行列の演算
行列には次のような演算が用意されています．
+ 行列の和 - 行列の差 determinant(A) 行列式 transpose(A) 転置行列 A^n 行列の成分ごとのベキ乗 A^^n 行列のベキ乗 (n=-1 で逆行列) A . x 行列 A とベクトル x の積 A . B 行列 A, B の積
行列 A の行列式，転置行列，逆行列をそれぞれ求めてみましょう．
(%i2) determinant(A); (%o2) - 2 (%i3) transpose(A); [ 1 - 3 ] (%o3) [ ] [ - 2 4 ] (%i4) A^^-1; [ - 2 - 1 ] [ ] (%o4) [ 3 1 ] [ - - - - ] [ 2 2 ]
18.3 行列と行列・ベクトルの積
行列 A と行列 B の積は A . B ，
行列 A とベクトル x の積は A . x で計算できます．
(%i5) A . A^^-1; [ 1 0 ] (%o5) [ ] [ 0 1 ] (%i6) x: [1, 2]; (%o6) [1, 2] (%i7) A . x; [ - 3 ] (%o7) [ ] [ 5 ]
18.4 問題
- 行列 A とその逆行列の積が単位行列になることを確かめよ．
- 連立方程式
x + y = 2,2*x - y = 1を逆行列を利用して解け． （ヒント：係数行列はmatrix([1, 1], [2, -1])^^-1， 右辺のベクトルは[1, 1]と表される）
19 いろいろな練習問題
- 多項式 P, Q, R を
P = x^2 + x - 3,Q = 2*x^2 - x + 4,R = -3*x^2 + 5とする．P - Q + Rを計算せよ。 - 方程式
x^2 = 2の解を(1, 2)の範囲で求めよ． - 方程式
sin(x) = x-1の解を(0, %pi)の範囲で求めよ． - 二次関数
y = x * exp(-x^2)が極大値とその時の x の値を求めよ． 1/(1+x^2)の不定積分を求めよ．1/(1+x^2)を[0, 1]の範囲で積分せよ．x*sqrt(x-1)の不定積分を求めよ． ヒント：ratsimp(),factor()で式を簡単化する．sin(x+%pi/4)を sin(x), cos(x) の和で表せ．sin(3*x)を sin(x) だけで表せ． ヒント： 式 eq 中の A を B で置き換えるにはratsubst(B, A, eq)とする．tan(x)をx=0の周りで5次の項までテイラー展開せよ．sin(x)/xをx=0の周りで5次の項までテイラー展開せよ．