« 整数演算のみで楕円を描く | メイン | Maximaで微分方程式 »

 Maximaでドラゴンカーブ

MathematicaとかMapleとかの数式処理ソフトは、何であんなに高いのだろう。2008/7/26現在、個人向けライセンス価格はMathematicaは40万円超え、Mapleも30万円超えだ。教育機関にも政府機関にも属していない、日曜数学ファンの私にはとても手が出せない。
会社でも研究職とかでなく、ソフトウェア開発の付帯業務をやってるので、数学は一切使っておらず、買ってもらえる見込みは無い。バグ収束曲線が云々で必要とか言っても、そんな統計処理ならExcelで十分と言われてしまうであろう。
従って、MathematicaやMapleに触れる機会が無い。

そんな折、Maximaというフリーの強力な数式処理ソフトがあることを知ったので、今回ヤボ用があって、インストールして使ってみた。
とりあえず、テーラー展開させてみる。

 KURO-BOX> maxima

Maxima 5.10.0 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.7 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
This is a development version of Maxima. The function bug_report()
provides bug reporting information.
(%i1) taylor(sin(x),x,0,10);
                          3    5      7       9
                         x    x      x       x
(%o1)/T/             x - -- + --- - ---- + ------ + . . .
                         6    120   5040   362880
(%i2) taylor(sqrt(r^2-y^2),y,0,10);
                    2      4      6         8       10
                   y      y      y       5 y     7 y
(%o2)/T/       r - --- - ---- - ----- - ------ - ------ + . . .
                   2 r      3       5        7        9
                         8 r    16 r    128 r    256 r
(%i3)
まあ言いたいことは十分わかるが、表示がなんか残念だ。
そこで、満足のいく表示を得るためにTeXmacsを使ってみる。以下の画面はMaxima 5.11.0とTeXmacs 1.0.6.11を使って実際に表示されたものである。%i??という青い部分が入力した部分だ。(表示環境はFreeBSD 6.3+Xorg 1.4.0)

積分もしてみる。

e^(-x^2)を-∞から+∞まで積分すると√πになるのは、手計算でやると面積分にして極座標変換してと技を尽くして偶然解ける感じの有名な例題だが、Maximaは区分求積による近似値計算ではなく、あっさり解析的に計算したようだ。

せっかくTeXmacsを使ってるので、見た目が派手な連分数を表示してみる。

cf(continued fraction)の出力は連分数のリスト表記になっていて、それをcfdisrepすると連分数表記になる。それをratで実数にすると(無限の連分数でないので元々実数だが)約分され、floatで小数にすると元の√2に限りなく近いことがわかる。
蛇足だが、cfdisrepしないと実数や小数に変換できない。

√3についてもやってみる。

方程式を解いてみる。solveで方程式を解析的に解ける。

解けない場合は、xの解にxが混ざったりする。(例は省略)

allrootsを使うと、多項式の方程式を数値的に計算できる。

指数の方程式を渡すと文句を言われる。

find_rootを使うと、(おそらくニュートン法による)近似計算で、解が1つだけ含まれる範囲の解を計算できる。

x^3=3^xのe~100の間の解は当然3なのだが、誤差が出た。


さて、せっかくフリーソフトなので、参考資料も無料のものを、と思って、図書館でMaximaの本を探したら、無かった。その代わり、あるMathematicaの本で、わずか十数行のコードでドラゴンカーブ(ドラゴン曲線)を描くMathematicaのコードを見つけた。マンデルブロー集合とかペアノ曲線とかのフラクタル図形の簡単なやつだ。Maximaの練習とMathematicaとの比較を兼ねて、その本を参考にMaximaで似たようなコードを書いてみた。

dragon(segments_):= block([p, q, r, s],
 p: map(
  lambda([x], [
   [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
   [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
  ]), segments_),
 q: map(
  lambda([x], [
   {[x[1][1], x[1][2]], [x[1][2], x[1][3]]},
   {[x[2][1], x[2][2]], [x[2][2], x[2][3]]}
  ]), p),
 r: map(lambda([x], listify(x)), flatten(q)),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(r)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 r
)$

p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$
for i:1 thru 10 do p:dragon(p)$
gnuplotが使える環境でこれを実行すると、以下のようなドラゴンカーブが描かれる。TeXmacsは必要ない。順に、L字型に配置した2本の線から始めて、全体を90°回転させて先端に繋げる、というのを1回、2回、…10回繰り返したものである。
アルゴリズムとしては、全ての線をL字に折り曲げる、というのを繰り返せばドラゴンカーブと同じ形になるというのを使っている。必要最小限にコードを簡略化して難解になったし、同様のアルゴリズムの説明は多くあると思うので、ここではコードの説明は省く。

段階が進むと上記コードの処理時間がやたら長くなるので、Perl的発想で、mapしまくって余計なリストを作ってるからかと思って、なるべくリストを作らないように、以下のような2つのコードを作って試してみた。

dragon2(segments_):= block([p, q, x, s],
 p: map(
  lambda([x], [
   [x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
   [x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
  ]), segments_),
 q: [],
 for x in p do q: append(q, [
   [[x[1][1], x[1][2]], [x[1][2], x[1][3]]],
   [[x[2][1], x[2][2]], [x[2][2], x[2][3]]]
  ]),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(q)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 q
)$
dragon3(segments_):= block([p, p1, p2, x, s],
 p: [],
 for x in segments_ do (
  p1: x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2,
  p2: x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2,
  p: append(p, [
   [[x[1][1], p1], [p1, x[1][2]]],
   [[x[2][1], p2], [p2, x[2][2]]]
  ])
 ),
 s: map(lambda([x], [realpart(x), imagpart(x)]), flatten(p)),
 plot2d([discrete, float(s)],
  [gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
 p
)$
結果は、少しずつは速くなったが、18.9秒が18.6秒とか18.0秒になるくらいで、ほとんど変わらなかった。余計なことはせず全て任せとけと言われた気分になった。Lispのリストの作成と削除は重くないらしい。CやPerlのような手続き型言語とは特性が違うのかも知れない。

さて、Mathematicaと比較すると、Maximaはわかりやすい参考資料が少ないのが非常につらい。taylorやintegrateの説明は当たりがつくのですぐに見つかったが、何かにつけ、どこを探せばいいのかがわからない。
今回MathematicaのNSolveに相当するfind_rootを見つけるのにもかなり時間がかかった。MathematicaのFlattenに相当するものが見つからず(Maximaのflattenは何重のリストでも全展開になる)、ドラゴンカーブのコードを実現するために数時間かけて解決策をひねり出したが、スマートなコードにならなかった。gnuplot関係はカットアンドトライでなんとか動いたが、結局わからないことだらけだ。
また、エラーメッセージが一昔前のCコンパイラ並に不親切で、エラーが出ると根本原因を見つけるのに一苦労でうんざりする。

しかし、MathematicaもMaximaも深くは知らないが、数式処理に関しては今の所MathematicaにできてMaximaにはどうしてもできないようなことは見当たらず、Maximaは非常に強力だと思った。Lispで書かれていることと合わせて、その出来栄えはフリーソフトとしては驚異的だと思う。

トラックバック

このエントリーのトラックバックURL:
http://ynomura.dip.jp/cgi-bin/mt/mt-tb.cgi/91

コメント (12)

GO_MAXIMA:

Barton のflatten.lisp を改良してMathematicaのように使えるようにしてみました。
こんな感じ
(%i4) load("newflatten.lisp");
(%o4) newflatten.lisp
(%i5) nflatten([[[a,b,c],d],[e,f]],1);
(%o5) [[a, b, c], d, e, f]
(%i6) nflatten([[[a,b,c],d],[e,f]]); /*default level 3*/
(%o6) [a, b, c, d, e, f]
(%i7) nflatten([[[a,b,c],d],[e,f]],2);
(%o7) [a, b, c, d, e, f]
関数定義はいかのとおりです。
(defun $nflatten (e &optional (n 5))
(setq e (ratdisrep e))
(cond ((or ($atom e) ($subvarp e)(or (member ($inpart e 0) (list '&^ '&=))))
e)
(t
(let ((op (multiple-value-list (get-op-and-arg e))))
(setq e (cadr op))
(setq op (car op))
(setq e (mapcar #'(lambda (x) (flatten-op x op n)) e))
(setq e (reduce #'append e))
(cond ((and (consp (car op)) (eq (caar op) 'mqapply))
(append op e))
(t
`(,op ,@e)))))))

(defun flatten-op (e op nlev)
(let ((e-op) (e-arg))
(setq e-op (multiple-value-list (get-op-and-arg e)))
(setq e-arg (cadr e-op))
(setq e-op (car e-op))
(cond ((and (>= nlev 1)(equal e-op op))
(mapcan #'(lambda (x) (flatten-op x op (1- nlev))) e-arg))
(t
(list e)))))


ynomura:

ありがとうございます。
最初のコードは、安全でない{}(set)とmap1つが無くなって、少しシンプルになりました。
dragon(segments_):= block([p, q, r],
p: map(
lambda([x], [
[x[1][1], x[1][1] + (x[1][2] - x[1][1]) * (1-%i)/2, x[1][2]],
[x[2][1], x[2][1] + (x[2][2] - x[2][1]) * (1+%i)/2, x[2][2]]
]), segments_),
q: nflatten(map(
lambda([x], [
[[x[1][1], x[1][2]], [x[1][2], x[1][3]]],
[[x[2][1], x[2][2]], [x[2][2], x[2][3]]]
]), p), 1),
r: map(lambda([x], [realpart(x), imagpart(x)]), nflatten(q,2)),
plot2d([discrete, float(r)],
[gnuplot_preamble, "set xrange [-1.2:1.2];set yrange [-1.2:1.2]"]),
q
)$
p: [[[0, 1/2-%i/2], [1/2-%i/2, 1]]]$
for i:1 thru 5 do p:dragon(p)$

なお、ご提示頂いたnflatten()を動かすには、flatten.lispにあるget-op-and-arg()を明示的に定義する必要があるようです。
当方の環境では、load(newflatten);の前にload(flatten);すると動きました。

gb:

連分数で漂着致しました。

http://f.hatena.ne.jp/mathnb/20100528021239
(内部リンク:20100528021239.gif

この 廣大の函数fが    ◆  いい加減法 (と命名します);

   x^2=7

3倍し;3x^2=3*7

  8*xを(いい加減)加え

3x^2+8*x=3*7+8*x

x*(3*x+8)=8*x+21

から 生まれた。なんて 信じる 学習者は 世界に 存在しない。

授業で いい加減法で 導出される方 は 存在しそう(嗚呼)......◆

★★ 廣大の函数f の導出過程を ご教示ください★★

(f の 導出にこそ 意味が在ると 考えます ので) 

---------------------------------------------------------------------

           また 
http://f.hatena.ne.jp/mathnb/20100528021239
    に倣い 例えば
Sqrt[3], Sqrt[109], Sqrt[263], Sqrt[431], Sqrt[601],
Sqrt[773], Sqrt[971], Sqrt[1153]
     等のそれぞれについて
廣大の函数f に相当する函数の導出を、 遊び心で、お願い致します;

f(Sqrt[3])=Sqrt[3](不動点) f[x]=
f(Sqrt[109])=Sqrt[109](不動点) f[x]=

ynomura:

連分数とは関係ないですね。
問題を解いてみればわかりますが、Nを問題で言う所の7として、f(Sqrt[N])=Sqrt[N]でf(x) > xでf'(x) > 0であれば何でもいいのです。(収束が速いか遅いかだけの違い、件の広大の入試問題からはそれ以上の条件は読み取れませんでした)
f(x)の形を(bx+c)/(ax+b)に限定し、a,b,cが全て整数、b^2 - a^2 * N = 1という制約を付けると、それは整数問題なので、総当たりで計算するしか無いと思います。

ちなみに、
N=3: a=1, b=2, c=3 (→ f(1)=1.667, f(f(1))=1.732)
N=263: a=8579, b=139128, c=2256277
N=109, 431, 601等については10^6までに適当な整数の組み合わせはありませんでした。
N=107ならa=93, b=962, c=9951
がありました。

b^2 - Na^2 = 1という制約を外すと、条件はc=aNとb^2>N*a^2だけになりますので、例えば広大の問題の8を9以上に変えても良いことになります。
N=109: a=3, b=32, c=109*3
N=431: a=3, b=63, c=431*3
N=601: a=3, b=74, c=601*3
(以下略)
いかがでしょうか。

gb:

同じ文面で たずねた ところ 或る方から
連分数
    [2,1,1,1,2+x]=2+1/(1+1/(1+(1/(1+1/(2+x)))))=(8*x + 21)/(3*x + 8)
が広島大学で出題されたf(x)の正体。
と 返事をいただきました。

連分数数と大いに関係が在ります。

Sqrt[61]等で f を 導出願います。

ynomura:

面白い情報ありがとうございます。
そういう連分数で書ける式に限定する条件、またはそれに等価な条件は、件の広大の問題には含まれていないと思いますが、まあ形からして元ネタはそれなんでしょうね。
Maximaでやってみました。

f(N) := block([y],
y:cf(sqrt(N)),
y[length(y)] : x+floor(sqrt(N)),
ratsimp(cfdisrep(y)));

f(7);
→(8*x+21)/(3*x+8)
f(61);
→(29718*x+232105)/(3805*x+29718)
f(109);
→(8890182*x+92816225)/(851525*x+8890182)

gb:

再読し 気付き  追記します;
http://www.wolframalpha.com/input/?i=e%5E%28-x%5E2%29

gb:

>Maximaでやってみました。
多謝!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

mathematica に ヤラセるには
如何に記述すればよいのでしょうか?
ご教示下さい。

ContinuedFraction[Sqrt[N] 云々....

gb:

再読し 追記;
http://www.wolframalpha.com/input/?i=z%5E17%3D1

0の原像 f^(-1)(0) のみでなく
fによる像 f[D] etc
mathematicaで叶いmath.

ynomura:

Mathematicaだと、こんな感じでヤラセられると思います。
F[n_] := Module[{y = Flatten[ContinuedFraction[Sqrt[n]]]},
y[[Length[y]]] = x + Floor[Sqrt[n]];
Simplify[FromContinuedFraction[y]]
]
F[7]
F[61]

いつの間にか、MathematicaのWebインターフェイスができてたんですね。
見よう見まねで、一応そのWolframAlphaでもやってみました。
http://www.wolframalpha.com/input/?i=Simplify%5BFromContinuedFraction%5BReplacePart%5BFlatten%5BContinuedFraction%5BSqrt%5B%23%5D%5D%5D%2C+-1+-%3E+x+%2B+Floor%5BSqrt%5B%23%5D%5D%5D%5D%26+%2F%40+%7B7%2C41%2C61%7D%5D
残念ながら1行1式でしか書けないようですので、逐次処理のプログラミングに慣れた者には結構厳しいですね。

gb:

In[62]:=
Table[Prime[k], {k, 4, 19}]

Out[62]=
{7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67}

In[64]:=
{F[Prime[4]], F[Prime[13]], F[Prime[18]]}

Out[64]=
{(21 + 8*x)/(8 + 3*x), (205 + 32*x)/(32 + 5*x), (232105 + 29718*x)/(29718 + 3805*x)}

WolframAlphaは更に或点で展開表示していますが「なんでだろー」

ynomura:

b^2 - a^2 * N = 1を満たす整数a,bを求める問題は、「ペル方程式」と呼ばれる、色々な所で出てくる方程式だそうです。上記の広大の問題は、やはり連分数と直接の関係は無く、「ペル方程式」を応用して作られた問題だと考えられます。
また、言うまでもなく、この問題はニュートン法の問題として作られたものであり、導出過程と出題意図とは関係ありません。

f(x)の形を(bx+c)/(ax+b)に限定し(平方根を分数式で近似することを考える、分子を(cx+d)としても結局b=cになる)、a,b,cが全て整数、f(√N)=√Nとすると、このfを繰り返し適用する時に収束が速くなるには、x=√Nの近傍でf'(x)=(b^2-aN)/(ax+b)^2の絶対値がなるべく小さい方が良い(大きいほどf(x)が√Nから遠ざかる)ので、b^2 - a^2 * N = ±1という制約がつきます。従って、広大の問題のa,bは「ペル方程式」の最小解のことだと思います。(連分数でも求まるが、連分数を使わなくても求まる)
f(x)を繰り返し適用すると、分母にペル方程式の解が次々に出てきますが、これもペル方程式の特性を利用したものと考えられます。

(参考)271828さんの滑り台Logへのリンク
ペル方程式と共役の考え
無理数と連分数(追記あり)
残念なことに、gbさんのコメントは、同一文面を多方面に送信していたマルチポストだったようです。何人かの方が個別に同じ問題に取り組んでおられました。無駄な労力を使わされたようで、あまり気分の良いものではありません。

コメント投稿フォーム

※投稿されたコメントはオーナーが承認するまで表示されません。


Powered by Movable Type 3.35