04 Maxima,Rと連携

Maximaと連携

KeTCindyでは,数式処理ソフトのMaximaと連携することができます。これにより,方程式を解いたり,微積分の計算を行ったりすることが数式処理としてできるようになります。

基本的な使い方

Maximaと連携をするための関数の基本的な使い方です。

Maxfun(name,式,リスト,option)

Maximaの関数を実行してその結果を取得します。第2引数の「式」がMaximaの関数名,第3引数の「リスト」が,その関数に渡す引数のリストです。

これを実行すると,バッチファイル(Macではシェルファイル)でMaximaが実行され,その結果を Fhead.txt に書き出し,これをKeTCindyが読み込みます。(Fheadは書き出すファイル名の頭部として設定されているものです)読みこんだ結果は,変数 mx+name に代入されるとともに,コンソールに表示されます。optionに "Disp=no"をつけると,コンソールでの表示はしません。

また,この関数の戻り値は,変数 mx+name に代入されたものと同じで,第1引数の式に1つでも文字があると文字列となり,すべて数字(+,-, .

を含む)の場合は16 桁以下であれば数,それ以上の場合は文字列となります。

例:関数の微分

Maximaで関数の微分を行うには,diff() を使って次のようにします。

diff(sin(x),x)

これで,sin(x)の導関数 cos(x)が求められるのですが,これを使ってたとえば微分係数を求めるには,次のようにする必要があります。

define(g(x),diff(sin(x),x));

すると,g(%pi/3) で cos(π/3) の値が求められます。

CindyscriptでもMaximaでも,関数の定義は g(x):=・・ のようにするのですが,ここでは 上のようにやらなければうまくいきません。

g(x):=diff(sin(x),x)

ではだめなのです。

これをKeTCindyでやるには次のようにします。

Mxfun("1", "diff",["sin(x)","x"]);

g(x):=parse(mx1);

で g(x) が cos(x)として定義されます。

1行で g(x):=parse(Mxfun("1", "diff",["sin(x)","x"]));

とすることもできます。

Maximaでは第n次導関数も簡単に求められます。

diff(sin(x),x,2)

とすれば第2次導関数が求められます。

KeTCindyでは

g(x):=parse(Mxfun("1", "diff",["sin(x)","x","2"]));

とすればよいですね。

Maximaでの数式処理と,KeTCindyでの書式をいくつか例示しておきましょう。

有理式(分数式)の通分 ratsimp(x/(x+2)+(x+3)/(x+1)); Mxfun("1","ratsimp",["x/(x+2)+(x+3)/(x+1)"])

通分した式の分子だけとり出す num(ratsimp(x/(x+2)+(x+3)/(x+1)));

Mxfun("2","num",["ratsimp(x/(x+2)+(x+3)/(x+1))"]);

通分した式の分母だけとり出す denom(ratsimp(x/(x+2)+(x+3)/(x+1)));

Mxfun("3","denom",["ratsimp(x/(x+2)+(x+3)/(x+1))"]);

因数分解 factor(x^3-1);

Mxfun("4","factor",["x^3-1"]);

式の展開 expand((x^2+x+1)*(x^2-x+1));

Mxfun("5","expand",["(x^2+x+1)*(x^2-x+1)"]);

整式A÷Bの商を求める quotient(x^3+x^2+2,x+3,x);

Mxfun("6","quotient",["x^3+x^2+2","x+3","x"]);

整式A÷Bの余りを求める remainder(x^3+x^2+2,x+3);

Mxfun("7","remainder",["x^3+x^2+2","x+3"]);

整式A÷Bの商と余りを求める divide(x^3+x^2+2,x+3);

Mxfun("8","divide",["x^3+x^2+2","x+3"]);

部分分数への展開 partfrac((x^2+2*x+3)/((x+3)*(x-1)),x);

Mxfun("9","partfrac",["(x^2+2*x+3)/((x+3)*(x-1))","x"]);

方程式を解く solve(x^3-2*x^2-x+2=0,x);

Mxfun("10","solve",["x^3-2*x^2-x+2=0","x"]);

※solve(x^3-2*x^2-x+2,x); としても,方程式と解釈して解いてくれます。

連立方程式を解く solve([y=x^2,y=x+2],[x,y]);

Mxfun("11","solve",["[y=x^2, y=x+2]","[x,y]"],[""]);

不定積分を求める integrate(x*sin(x),x);

Mxfun("12","integrate",["x*sin(x)","x"]);

CalcbyM(name,リスト)

Maximaのコマンドを実行します。Mxfun()との違いは,1つのコマンドではなく,複数のコマンドからなる一連の処理を実行できることです。リストにその一連の処理を,関数,引数 の2つずつをセットにしてそれぞれ文字列にして記述します。

例:関数の微分係数を求める

Maximaで関数の微分を行ない,微分係数を得るには次のようにするのでした。

define(g(x),diff(sin(x),x));

g(%pi/3);

ここで,g(x) を導関数とするかわりに,式の評価をする ev() という関数を使って微分係数を求めることができます。

gx:diff(sin(x),x);

ev(gx,x=%pi/3);

とすると,gx に diff(sin(x),x) の結果,すなわち cos(x) が代入され( : が代入文です) ev(gx,x=%pi/3); によって x を π/3 として「評価」します。つまり,cos(x)のxに %pi/3 を代入したのと同じ結果を得るのです。

これをKeTCindyで行うには

cmdL=[

"gx:diff",["sin(x)","x"],

"c:ev",["gx","x=%pi/3"],

"c",[]

];

CalcbyM("tn1",cmdL);

とします。

ここでいくつか注意があります。

まず,Maximaでの実行と同様な "ev",["gx","x=%pi/3"] では結果が得られないことです。"c:ev" として,変数 c に代入する必要があります。

次に,戻り値を得るために,最後にその変数と空リストをセットにして与えます."c",[] がそうです。すると,cの値が name の tn1 に代入されます。

Mxfun()は戻り値と変数の値が同じでしたが,CalcbyMでは関数の戻り値はありません。(未定義)

したがって,上記のようにして戻り値を取得します。

Mxtex(name,式)

式をTeX 書式にします。

第2引数の式は,直接書いた式もしくはMxfun の戻り値。これをTeX の書式にします。

戻り値は,変数txname にも代入されm。

例:部分分数への分解

Mxfun("1","partfrac",["x^3/((x+1)*(x+2))","x"]);

Mxtex("1",mx1);

Expr([0,1],"e",tx1);

ここで,mx1,tx1 はそれぞれMxfun("1",・・) , Mxtex"1",・・) の結果(戻り値)です。

mx1,tx1 はコンソールにも表示され,tx1 は次のようになっています。

\frac{8}{x+2}-\frac{1}{x+1}+x-3

Cindyscript はTeX 書式をサポートしているのでこれで描画面に数式が表示されます。そして,それはそのままTeXのテキストに出力できます。

2次関数のグラフとx軸との交点の座標

Maximaとの連携で,Cindyscriptだけでは面倒,もしくはできない処理ができるようになります。そんな利用例です。

2次関数のグラフを描き,x軸との交点の座標を表示します。

Cinderellaを使う利点は,インタラクティブに図形を動かせることです。2次関数のグラフも,頂点をドラッグしてグラフの位置を変え,方程式を表示することができます。ここまではCindyscriptで簡単にできますが,このときx軸との交点の座標を表示しようとするとうまくいかないことがあります。頂点がどこにあっても2次方程式の解の形で表示するのに,Maximaによる2次方程式の解の取得が役に立ちます。

まず,Cinderellaの作図ツールで点Aをとっておきます。これが放物線の頂点です。

スクリプトは以下の通りです。

=============================

if(A.y<0,

fx="(x-"+text(A.x)+")^2"+guess(A.y),

fx="(x-"+text(A.x)+")^2+"+guess(A.y);

);

drawtext([3,2],"y=$"+fx+"$");

cmdL=[

"ans:solve",[fx,"x"],

"ans",[]

];

CalcbyM("ans",cmdL);

p1=indexof(ans,"[");

p2=indexof(ans,",");

p3=indexof(ans,"]");

s1=substring(ans,p1,p2-1);

s2=substring(ans,p2,p3-1);

s1=replace(s1,"x =","");

s2=replace(s2,"x =","");

Mxtex("1",s1);

Mxtex("2",s2);

Plotdata("1",fx,"x");

Expr([parse(s1),-0.5],"e",tx1);

Expr([parse(s2),-0.5],"e",tx2);

=============================

1〜4行目で,点Aが頂点となる関数の式を作ります。

5行目で式を画面に表示します。(TeXのテキストには出力されません)

6行目でコマンドを定義し,CalcbyM("ans",cmdL) によって,2次方程式の解がリスト ans に出力されます。

このリストは文字列になっているので,そこから,解の部分をとり出し(s1,s2)

さらに,解の部分から x = を削除して,あらためて s1 s2 とします。

このとき,s1は -(sqrt(6)-2)/2 のような文字列です。

これを Mxtex で \frac{2-\sqrt{6}}{2} に変換します。TeX書式です。

あとは,グラフと解の表示です。

点Aをドラッグすると,インタラクティブにグラフと交点の座標が変化します。

ただ,Drawスロットに書いているので,ちょっと点を動かしただけでバッチ/シェルが実行されるので,リアルタイムではありますが反応が遅くなります。

これらを関数化してInitialization スロットに置き、ボタンを用意すれば,Aの位置を決めた後ボタンで実行するようにすることもできます。

なお,2次方程式はCindyscriptでも解けますが,数値計算で解きますのでそのままでは解の公式の形にはなりません。Maximaを使う方が便利でしょう。

テイラー展開

関数のテイラー展開をして,もとの関数のグラフと比べます。

Cindyscriptでも微分はできますが,数値計算をしているので4次くらいまでが限界で,5次になると誤差が大きくなってしまいます。Maximaなら何次でもできますし,それも引数に次数を与えるだけでできます。(Cindyscriptでは微分の繰り返しにしなければなりません)

次は sin x のテイラー展開の例です。

Mxfun("1","taylor",["sin(x)","x",0,7],[""]);

Plotdata("1","sin(x)","x",["da"]);

Plotdata("2",mx1,"x");

Expr([[1,2],"e",Mxtex("1",mx1),[pi,-0.5],"s","\pi",[-pi,-0.5],"s","-\pi"]);

微分の次数をスライダで決めるようにすると,テイラー展開した関数のグラフが,sin x のグラフに近づいていく様子をインタラクティブに観察できます。

スライダの左端をy軸上に置き,1行目を次のように書き換えます。

n=ceil(C.x*2);

drawtext(C.xy-[0.3,0.5],"n="+n,size->20);

Mxfun("1","taylor",["sin(x)","x",0,n],[""]);

Cindyscriptだけではここまではできません。Maximaを利用する価値のある例です。

Rと連携

KeTCindyでは,Rと連携して統計計算を行うことができます。

これにより,箱ひげ図を描いたり,散布図を描いたりすることが簡単にできるようになります。

Boxplot(名前, データ, 垂直位置, 高さ,option)

箱ひげ図を描きます。データは変数に代入して与えます。

外部ファイルのデータを読み込むには,Readcsv()を使います。

例:表示領域の対角点A,Bを作図しておきます。(Aが左下)

Readcsv("1","datafile.csv");

dt1=apply(rc1,#_1);

dt2=apply(rc1,#_2);

Boxplot("1",dt1/20,1,1/2);

Boxplot("2",dt2/20,3,1/2);

Framedata2("1",[A,B]);

Rulerscale(A,["r",0,6,1],["f",1,"\mbox{dt1}",3,"\mbox{dt2}"]);

CalcbyR(変数名, コマンド列,option)

Rのコマンド列を与えて実行します。

Histplot(name,data,option)

ヒストグラムを描きます。

PlotdataR(name, 式, 変数)

R の関数のグラフを描きます。Cindyscriptの組込関数にないものを描けます。

PlotdiscR(name, 式, 変数)

dbinom (二項分布),dpois(ポアソン分布),dgeom(幾何分布)など離散型確率分布のグラフを描きます。

Scatterplot(name,lename,option)

2次元データを読み込み,散布図を描きます。回帰直線を描くこともできます。