pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

入力例で学ぶMaximaの使い方(入門)

数式変形をすることそれ自体も楽しいものですが、単純作業という一面があることも否めません。
効率よく勉強するには面倒な計算はコンピュータに任せてしまうのも一つの方法でしょう。 

数式を処理するソフトウェアとしてよく知られているのはMathematicaでしょうか?
利用者も多いので、わからないことがあったときなどは、その充実したヘルプ機能に加えて、
ネット検索によっても解決策が得られやすいという長所があります
しかしいかんせん有料(結構高い)なので個人利用で導入するのは気乗りしません。

ここではMathematicaの代替として遜色ないフリーソフトである『Maxima』の使い方を紹介します。
といっても、いち利用者として具体例を扱いながら、習うより慣れろで使い方を学んだだけなので、
ちょっと高度で便利な関数電卓の(自己流の)使い方を紹介します、という感覚の記事です。

Maximaのダウンロードやインストールについては、以下のサイト

maxima.sourceforge.net

や『maxima ダウンロード』などとネット検索して出てきたサイトを参考にしてみてください。



Maximaの使い方入門

Maximaの入力はセミコロン" ; "、またはドル記号" $ "で終わります。
後者の" $ "は結果を画面に出力させないためのオプションです。
評価したい式を入力したら『Shift+Enter』で計算を実行してくれます。

何はともあれ、まずは実際に使ってみましょう
以下では具体的なMaximaの入力例を示します。
コピー&ペーストして計算を実行させたり、入力をいじったりして遊んでみてください。

加減乗除・冪乗・階乗・二重階乗・ルート

3+2;   /* コメント(評価の際に無視される文字列)はスラッシュ記号とアスタリスク記号で囲みます */
3-2;
3*2;
3/2;
3^2;
50!;
50!!;
sqrt(50);

浮動小数点数

大きな数
float(10^8);
bfloat(10^8);
Napier数(自然対数の底
%e;   /* Napier数は" e "にパーセント記号をつけます */
float(%e);
bfloat(%e);
円周率
%pi;   /* 円周率は" pi "にパーセント記号をつけます */
float(%pi);
bfloat(%pi);
無限大
inf;
ルート2
float(sqrt(2));
bfloat(sqrt(2));

最大値・最小値

max(sqrt(2),2,%pi);
min(sqrt(2),2,%pi);

整数部分を取り出す

round(sqrt(2));
round(-9/7);

Gauss記号

ある数に対して、それ以下の最大の整数値を返します。

fix(3);
fix(5/2);
fix(%pi);

四捨五入

整数部分の取り出し関数roundと浮動小数点数の表示する関数floatを組み合わせて

小数を好きな桁数のところで四捨五入して近似した数を得ることができます。

float(9/7);

とすると

 \quad \frac{9}{7}=1.285714285714286

ですが

float(1/10^8*round(10^8*9/7));

とすると小数点以下第8位のところで四捨五入して近似した

 \quad \frac{9}{7}\sim 1.28571429

を出力させることができます。

こういった操作は、Maxima数値計算する場合に使えます。

複素数

虚数単位
%i;   /* 虚数単位は" i "にパーセント記号をつけます */
%i^2;
複素共役
conjugate(2+sqrt(3)*%i);
conjugate(a+b*%i);   /* 何も指定しなければ文字は実数として処理されます */
実部・虚部
realpart(2+%i*sqrt(3));
imagpart(2+%i*sqrt(3));
絶対値
abs(2+%i*sqrt(3));
極表示
polarform(2+%i*sqrt(3));
分母の有理化
rectform(1/(2+%i*sqrt(3)));
rectform(1/(a+%i*b));

代数的操作

素因数分解する
factor(10!);
ifactors(10!);
整数の商と余りを求める
divide(16,3);

商だけ、余りだけを求めるには、それぞれ

quotient(16,3);
mod(16,3);
自然数を法とする合同
is(mod(6, 5)=mod(31, 5));
最大公約数
gcd(6,18);
最小公倍数
lcm(6,18);
因数分解する
factor(x^2-4*x+3);
factor(x^6-1);
Gauss整数係数多項式因数分解する
gfactor(x^2+1);

Gauss整数:実部と虚部がともに整数である複素数

部分分数分解する
partfrac(1/(x^2-4*x+3), x);
partfrac(1/(x^6-1), x);
通分する
xthru(1/(x-1)^2+1/((x-1)*(x-2)));
多項式を展開する
expand((x+1)^6);
式を簡単化する
ratsimp((2*x-2)/((x-1)^2*(x-2)));

ただし根号をともなう場合にはratsimpよりradcanのほうが有用かもしれません:

ratsimp((x-1)/sqrt(x^2-1));
radcan((x-1)/sqrt(x^2-1));

ちなみに、多項式の場合にはratsimpもradcanも同じ結果を与えます:

ratsimp((2*x-2)/((x-1)^2*(x-2)));
radcan((2*x-2)/((x-1)^2*(x-2)));

ですので、多くの場合ではradcanでじゅうぶんかもしれません。

違いが出るのは対数関数などをともなう場合です:

ratsimp(log(x/y));
radcan(log(x/y));
式の真偽を判定する
is(0=sqrt(2)-2^(1/2));
is(0=sqrt(2)-2^(1/3));

ただし数式を適当に処理しないとうまく判定しない場合も多いです: 

is(x^2-2*x+1=(x-1)^2);
is(x^2-2*x+1=expand((x-1)^2));

また、等しくないことを確認する場合には" = "の代わりに" # "を使います:

is(1#1);
is(1#0);
多項式から特定の項の係数を取り出す
ratcoeff((x + y)^6, x, 3);
ratcoeff((x + y)^6, x, 4);
値を代入する
subst(3, x, a*x^2+b*x+c);
総和・総乗
sum(k, k, 1, 10);
sum(k, k, 1, n), simpsum;
prod(k, k, 1, 10);
is(10!=prod(k, k, 1, 10));

方程式

連立方程式を解く
linsolve ([3*x+4*y=5, 2*x+3*y=3], [x, y]);
代数方程式の根を求める
solve(a*x^2+b*x+c=0, x);

以下の記事ではMaximaを用いた代数・連立方程式の計算についてまとめています:

pianofisica.hatenablog.com



極限

極限値を求める
limit( sin(x)/x, x, 0 );

関数の定義・変数への代入

文字に式を割り当てる、値を代入する
u:x+y;
v:(x+y)^3;
subst(5, x, u);
ev(u, x: 5);
ev(v, [x: 3, y: 2]);
割り当てた文字を使う・四則演算
u+v;
expand(u+v);
factor(u+v);
expand(u*v);
割り当てた文字を使う・微分する
diff(u, x);     /* u の x についての1階偏微分係数 */
diff(u, y);     /* u の y についての1階偏微分係数 */
diff(v, y, 2);     /* v の y についての2階偏微分係数 */
is(0=diff(u, x, 1)-diff(u, x));
割り当てた文字を使う・積分する
integrate(u, x);   /* u を x について不定積分します */
integrate(u, x, x1, x2);   /* u を x について x1 から x2 までの区間で定積分します */
integrate(v, y);
関数を定義/微分積分/Taylor展開する・代入する
f(x):=log(1+x)/x$
integrate(f(x), x);
diff(f(x), x, 1);
taylor(f(x), x, 0, 5);  /* 関数 f(x) を x について x=0 の周りで5次までTaylor展開します */
g(x):= x^2$
h(x):=g(f(x))$
subst(5, x, h(x));
等式の右辺・左辺、式の分母・分子を取り出す
F(x):=x+1$
G(x):=x^2+x-2$
eq: x+1 = F(x)/G(x)$
rhs(eq);     /* eq の右辺 */
lhs(eq);     /* eq の左辺 */
denom(rhs(eq));     /* eq の右辺の分母 */
num(rhs(eq));     /* eq の右辺の分子 */
solve(eq, x);

数値積分

Maximaには数値積分の組み込み関数としてRomberg積分が用意されています。

romberg(sqrt(1-x^2), x, 0, 1);  /*romberg(expr, var, min, max)*/

これは半径1の円の第1象限部分の面積なので、円周率の4分の1

bfloat(%pi/4);

に等しいはずです。しかしながら、Romberg積分の精度はそれほど高いわけではなく、

小数点以下5桁めから値のズレが認められます。

微分方程式

常微分方程式を解く
assume(k>0)$
desolve( diff( u(t), t, 2 ) = - k^2*u(t), u(t)  );

以下の記事ではMaximaを用いた常微分方程式の計算についてまとめています:

pianofisica.hatenablog.com

三角関数双曲線関数)の諸操作

引数の和を展開する
trigexpand(sin(a+b));
trigexpand(cos(a+b));
trigexpand(cosh(a+b));
trigexpand(tanh(a+b));
簡単化する
trigsimp(sin(a)^2+cos(a)^2);
trigsimp(sinh(a)^2-cosh(a)^2);
簡約する
trigreduce(sin(a)^2);
trigreduce(cosh(a)^3);
指数関数で表示する
exponentialize: true$  /* 元に戻すには true を false に変えます */
sin(a);
cosh(a);

フーリエ級数展開する

詳しくは以下の記事を参照してください。
pianofisica.hatenablog.com

リスト

リストを作る・要素を参照する
S:[a, a, a, b, c];
S[1];
S[4];
リストを利用する
T:[1, 10, 100, 1000]$
sum( T[i], i, 1, 3);
複数のリストを結合する
U:append(S,T);
リストの先頭に要素を加える
cons(h, U);
リストの末尾に要素を加える
endcons(e, U);
要素の最大値・最小値
lmax(T);
lmin(T);
要素の個数
s:length(S);
t:length(T);
is(s+t=length(U));
リストを生成する
n[k] := k$
n[10];
N[n]:=makelist(n[k], k, 1, n)$
N[10];
集合を入力する
S0: {a, a, a, b, c};    /*  集合は波括弧で挟んで入力します */

上で入力したリストSとの違いを比較してください。集合は同一要素の重複を省きます。また、上記のリストへの操作、append、cons、endcons、lmax、lmin、lengthは集合に対しても同様に使えます。

和集合・共通部分をとる
S1: {a,b,c,d,e}$
S2: {a,b,f,g}$
union(S1,S2);
intersection(S1,S2);
集合をリストに変換する
listify(S0);
リストを集合に変換する
setify(S);



繰り返し操作

WHILE
i:0$ while i<5 do(i:i+1,print("i=", i));
FOR
for j:1 step 1 thru 5 do print(j);
応用:逐次代入による方程式の数値解

方程式  x+e^x=0  x_0=0 として逐次的に数列

 \displaystyle{\quad x_{n+1}=-e^{x_n} }

に代入することで、数値的な近似解を求めることができます:

kill(all)$
x[0]:0$
for i:1 step 1 thru 100 do x[i]:bfloat(-exp(x[i-1]))$
x[100];
bfloat(x[100]+%e^(x[100]));

あるいは while 文を用いると

kill(all)$
x[0]:0$
i:0$ while i < 100 do(
i:i+1,
x[i]:bfloat(-exp(x[i-1]))
)$
x[100];
bfloat(x[100]+%e^(x[100]));

TeX形式で計算結果を出力する

Eq:a*x^2+b*x+c=0$
sol:solve(Eq, x)$
sol1:rhs(sol[1]);
sol2:rhs(sol[2]);
tex(sol2);

計算にかかった時間を表示する

showtime:true$

多項式の昇ベキ・降ベキ表示を切り替える

powerdisp:true$
expand((x+1)^5);
powerdisp:false$
expand((x+1)^5);

条件分岐

H(x):=if  x < 0  then 0 else 1$

引数が非負のとき1を、負のときにはゼロを返す関数です。

これはHeavisideの階段関数と呼ばれます。

応用:二分法による方程式の数値解

裳華房数値計算』(柳田・中木・三村)を参考にして

二分法によって方程式  \cos(x)=x の数値解を求めます。

if-then-else構文、複数の処理を1つの処理とみなすblock構文の具体例です:

f(x):=bfloat( cos(x)-x )$
a[0]:bfloat(0)$
b[0]:bfloat(1)$
i:0$ while i<10 do(
if  f((a[i]+b[i])/2) < 0  then block(a[i+1]:bfloat(a[i]), b[i+1]:bfloat((a[i]+b[i])/2))
else block(a[i+1]:bfloat((a[i]+b[i])/2), b[i+1]:bfloat(b[i])),
i:i+1);
c[n]:=bfloat((a[n]+b[n])/2)$
makelist(a[i], i, 0, 9);
makelist(b[i], i, 0, 9);
makelist(c[i], i, 0, 9);
makelist(f(c[i]), i, 0, 9);
応用:多項式の次数を出力する
kill(all)$
Ord(P, x):=
block(
n:1, while n<31 do(
c: ratcoeff(P, x, 30-n),
k[n]: if c = 0 then -1 else 30-n,
n:n+1),
lmax(makelist(k[n],n, 1,30)) 
)$
Poly: (t+1)^5$ Ord(Poly, t);

この関数を使ったさらなる応用として以下の記事を紹介しておきます:

pianofisica.hatenablog.com



グラフの図示

2次元プロット
plot2d (x^3+x^2-3*x+1, [x, -2, 2])$
片対数目盛
plot2d (exp(3*s), [s, -2, 2], logy)$
等高線
contour_plot (x^2 + y^2, [x, -5, 5], [y, -5, 5])$

行列

行列を入力する
s1: matrix([0,1], [1,0]);
s2: matrix([0,-%i], [%i,0]);
s3: matrix([1,0], [0,-1]);
M: matrix([ a11, a12, a13 ], [ a21, a22, a23 ]);
行列の要素を参照する
M[1][3];

あるいは

M[1,3];
行列の和・積をとる
s1+s3;
s1.s3;
行列式を計算する
determinant(s2);
逆行列を求める
invert(s2);
転置行列を求める
transpose(s2);
行列の固有値を求める
eigenvalues(s2);
行列の固有値固有ベクトルを求める
eigenvectors(s2);


以下の記事ではMaximaを用いた行列の計算についてまとめています:

pianofisica.hatenablog.com

pianofisica.hatenablog.com

pianofisica.hatenablog.com

直前の出力結果を引用する

A: p*x^2+q*x+r$
subst(3, x, %);   /* 直前の出力はパーセント記号で次のコマンドへ引き継ぐことができます */

文字を初期化する

A: p*x^2+q*x+r;
kill(A);
A;

すべての入力を初期化する

kill(all)$


また関連記事では、Maximaの使い方の応用編として、数学や物理学のさまざまな具体的な問題の計算にMaximaを使っています。Maximaの入力も載せているので、参考にしてみてください。


キーワードMaxima