フラクタルの最近の記事

今度は、解 ω と ω2 を実軸に近付けて、3 つの解がなす三角形を頂角が 60 度より小さい二等辺三角形にしてみます。

60minus-a800-480x270.jpg-0.5±0.8i

60minus-a700-480x270.jpg-0.5±0.7i

60minus-a600-480x270.jpg-0.5±0.6i

60minus-a510-480x270.jpg-0.5±0.51i

60minus-a500-480x270.jpg-0.5±0.5i

60minus-a400-480x270.jpg-0.5±0.4i

60minus-a330-480x270.jpg-0.5±0.33i

60minus-a300-480x270.jpg-0.5±0.3i

異変は 0.1i 刻みでは見逃してしまいます。0.01i 刻みで見ると、-0.5±0.51i と -0.5±0.33i のとき原点付近に黒い穴が空くことがわかります (0.001i 刻みにするともっと多くの場所で穴が空きます)。-0.5±0.33i の前後を 10 倍に拡大して見てみましょう。

60minus-b338-480x270.jpg-0.5±0.338i

60minus-b337-480x270.jpg-0.5±0.337i

60minus-b336-480x270.jpg-0.5±0.336i

60minus-b335-480x270.jpg-0.5±0.335i

60minus-b334-480x270.jpg-0.5±0.334i

60minus-b333-480x270.jpg-0.5±0.333i

60minus-b332-480x270.jpg-0.5±0.332i

60minus-b331-480x270.jpg-0.5±0.331i

60minus-b330-480x270.jpg-0.5±0.33i

60minus-b329-480x270.jpg-0.5±0.329i

60minus-b328-480x270.jpg-0.5±0.328i

60minus-b327-480x270.jpg-0.5±0.327i

60minus-b326-480x270.jpg-0.5±0.326i

60minus-b325-480x270.jpg-0.5±0.325i

60minus-b324-480x270.jpg-0.5±0.324i

白は反復回数は大きいけれど収束した初期値、黒は収束しなかった初期値です。色は白から黒に、形は丸から団子に、ドラスティックに変化しています。

ニュートン法とフラクタル に書いた 3 次方程式

x3-1 = 0

の解は

  1. 1
  2. ω = -1+sqrt(3)·i2
  3. ω2 = -1-sqrt(3)·i2

の 3 つです。これらは単位円に内接する正三角形の頂点の位置にあります。

ここで、解 1 を原点に近付けて、3 つの解がなす三角形を頂角が 60 度より大きい二等辺三角形にしてみます。

60plus-99-480x270.jpg0.99

60plus-97-480x270.jpg0.97

60plus-90-480x270.jpg0.9

60plus-70-480x270.jpg0.7

どんぐりがほつれてしまいました。

ニュートン法とアフィン変換

| コメント(0)

方程式

f(x) = 0

に対するニュートン法の漸化式は

xn+1 = xn - f(xn)f'(xn)

です。ここで f(x) の x をアフィン変換 (拡大、縮小、回転、平行移動) した関数

g(x) = f(p·x+q)

を考えると、方程式

g(x) = 0

に対するニュートン法の漸化式は

xn+1 = xn - g(xn)g'(xn)

xn+1 = xn - f(p·xn+q)p·f'(p·xn+q)

p·xn+1+q = p·xn+q - f(p·xn+q)f'(p·xn+q)

となります。これは f(x) の x をアフィン変換するとニュートン法の xn の軌道も同じようにアフィン変換されることを意味しています。xn の軌道には初期値 x0 も含まれますから、2 つの方程式の解の位置を頂点とする多角形が相似ならば、その 2 つをニュートン法で解いてどの解に収束するかで初期値を塗り分けたときに現れる図形も相似になります。ニュートン法とフラクタル で書いた

f(x) = x3-1

の場合は単位円に内接する正三角形の頂点の位置に解があって、初期値を塗り分けるとどんぐりが連なったような模様が現れました。ということは、3 つの解を結ぶと正三角形ができる 3 次方程式ならば、正三角形が複素平面上のどこにあってどちらを向いていても、ニュートン法で解いてどの解に収束するかで初期値を塗り分けると必ずどんぐりが連なったような模様が現れるのです。例えば 3 つの解を 1+i, -1+isqrt(3), 1-isqrt(3) に配置すると、こうなります。

newtons-method-2-480x270.jpg

ニュートン法とフラクタル

| コメント(0)

微分可能な関数 f(x) で表される方程式

f(x) = 0

が与えられたとき、初期値 xn を決めて漸化式

xn+1 = xn - f(xn)f'(xn)

の計算を繰り返すことで解を求める方法がニュートン法です。

3 次方程式

x3-1 = 0

の場合、複素平面上の初期値 x0 から出発して漸化式

xn+1 = xn - xn3-13xn2

の計算を繰り返すと、xn

  1. 1 に収束する。
  2. ω = -1+sqrt(3)·i2 に収束する。
  3. ω2 = -1-sqrt(3)·i2 に収束する。
  4. どの解にも収束せず複素平面上を彷徨い続ける。
  5. ゼロ除算で計算できなくなる。

のいずれかになります。この結果によって複素平面上の初期値 x0 を塗り分けます。1 に収束した初期値を赤、ω に収束した初期値を緑、ω2 に収束した初期値を青、それ以外を黒で塗ると、こうなります。

newtons-method-1-480x270.jpg

画像の範囲は左下が -1.6-0.9i、右上が 1.6+0.9i、中央が 0 で、解の位置は右側の赤が暗くなっている辺りが 1、左上の緑が暗くなっている辺りが ω、左下の青が暗くなっている辺りが ω2 です。

もしも、常に初期値から最も近い解に収束するのであれば、赤と緑、緑と青、青と赤の境界線はすべて直線になるはずです。しかし、実際の境界線は直線でないどころか、そもそも境界線と呼べるものは存在しません。なぜなら、このどんぐりが連なったような模様はフラクタル (自己相似的) で、どこまで拡大しても赤と緑の間には青が、緑と青の間には赤が、青と赤の間には緑が挟まっていて終わりがないからです。これは「平面が複数の色で塗り分けられているのに色の境界線が存在しない」という奇妙な図形なのです。

FC2 ジャンパピンを抜いた 060turbo でも VRAM は MMU でスーパーバイザプロテクトされていてそのままではユーザモードでアクセスすることはできません。しかし、スーパーバイザプロテクトを外してやればユーザモードで VRAM にアクセスできるようになります。具体的には VRAM に対応する論理ページのページデスクリプタのビット 7 をクリアすればよいのです。この方法は MMU を搭載していて FC2 がカットされているマシンでなければ使えませんし、そもそも Human68k 上で動作するプログラムは簡単にスーパーバイザモードに移行できるのであまり実用的とは言えません。強いて言えばシステム領域を保護したまま VRAM だけユーザモードでアクセスできるようになるので VRAM を直接叩くプログラムのデバッグが楽になるかも知れません。

fc2test.c
/*
  fc2test.c
    ユーザモードでVRAMにアクセスしてみるテスト
    2007-11-23 M.Kamada
  動作環境
    FC2ジャンパピンを抜いた060turbo
    FC2がカットされていてSYS_STATの$F002が実装されていれば
    MMU付きのMC68030に換装したX68030や040turboでも動くかも
  コンパイル
    コンパイラはgcc真里子版またはgcc2を使用
    ライブラリはXCまたはlibcを使用
    > gcc -Wall -m68040 -O6 -o fc2test fc2test.c -liocs
  実行
    > fc2test
*/

#include <stdio.h>
#include <iocslib.h>

#if 1
#define R0 (-0.75)
#define I0 (0.0)
#define SCALE (1.25 / 256)
#define DEPTH 64
#else
#define R0 (-0.738)
#define I0 (0.188)
#define SCALE (0.01 / 256)
#define DEPTH 512
#endif

#define PAGE_OFFSET_SIZE 0x00002000
#define PD_S 0x00000080

/* SYS_STAT $F002 */
long sys_stat_f002 (void *address, long attribute) {
  register void *a1 asm ("a1") = address;
  register long d2 asm ("d2") = attribute;
  register long d0 asm ("d0");
  asm volatile ("move.w #$f002,d1\n\
        moveq.l #$ac,d0\n\
        trap #15" : "=d"(d0) : "a"(a1), "d"(d2) : "d1");
  return d0;
}

/* fromからtoまでのスーパーバイザプロテクトを解除・設定する */
void supervisor_protect (void *from, void *to, int mode) {
  long s = mode ? PD_S : 0;
  void *address;
  long attribute;
  for (address = from; address <= to; address += PAGE_OFFSET_SIZE) {
    attribute = sys_stat_f002 (address, -1);
    if (attribute != -1 && (attribute & PD_S) != s) {
      sys_stat_f002 (address, (attribute & ~PD_S) | s);
    }
  }
}

int main (int argc, char *argv[]) {

  unsigned short *vram = (unsigned short *) 0x00c00000;
  int i, n, x, y;
  double ci, cr, ti, tr, zi, zr;

  /* 768×512ドット16色モードにする */
  if (CRTMOD (-1) != 16) {
    CRTMOD (16);
  }
  G_CLR_ON ();

  /* パレットを黄色のグラデーションにする */
  for (i = 0; i < 16; i++) {
    GPALET (i, (i << 12) | (i << 7) | (i << 1));
  }

  /* VRAMのスーパーバイザプロテクトを解除する */
  supervisor_protect (vram, vram + 1024 * 1024 - 1, 0);

  /* マンデルブロ集合を描画する */
  for (y = 0; y < 512; y++) {
    ci = I0 + (256 - y) * SCALE;
    for (x = 0; x < 768; x++) {
      cr = R0 + (x - 384) * SCALE;
      zr = cr;
      zi = ci;
      for (n = 0; n < DEPTH; n++) {
        tr = zr * zr;
        ti = zi * zi;
        if (tr + ti > 4.0) {
          break;
        }
        zi *= zr;
        zi += zi;
        zi += ci;
        zr = tr - ti + cr;
      }
      vram[x + 1024 * y] = (n / (DEPTH >> 4)) & 15;
    }
    if (B_KEYSNS ()) {
      break;
    }
  }
  while (B_KEYSNS ()) {
    B_KEYINP ();
  }

  /* VRAMのスーパーバイザプロテクトを設定する */
  supervisor_protect (vram, vram + 1024 * 1024 - 1, 1);

  return 0;
}

「かぐや」関連

| コメント(0)

「かぐや」が搭載している空間分解能 10m の地形カメラと空間分解能 20m のマルチバンドイメージャがとらえた初画像が公開されました。月面は火星ほど変化に富んでいるわけではありませんが、大気がないので小さなクレーターもことごとく保存されていてまるでフラクタル画像を見ているようです。詳細な地形データは月の歴史を紐解く手掛かりとなるだけでなく、将来の月面基地建設のための資料としても役立つことでしょう。極地方の日の当たらないクレーターの底がどうなっているのか気になります。

060 on 040turbo

| コメント(0)

ラキッ!さん から。X68030 のハードディスクから 060 on 040turbo のベンチマークの結果を拾ってきました。

- FSHARPU  MES( 8):◆ 会議《 X68 ハードウェア  》 96/03/26 -
01028/01035 LDN01325  かまだ           RE:060 on 040turbo
( 8)   96/03/26 07:31  01026へのコメント

060試作基板に関するベンチマーク報告(浮動小数点演算)


○ベンチマーク内容

        内容            マンデルブロ集合の描画
        範囲            -1.2504+0.0080i~-1.2409+0.0085i
        分解能          512×512ピクセル
        最大ループ回数  20000回/ピクセル


○動作環境

●X68000 PRO
システムクロック        10MHz
FPU                     純正コプロセッサボード
プログラム              コプロセッサボード直接アクセス
                        メインループを全てアセンブラで記述,16回分展開

●040turbo
システムクロック        25MHz
メモリモード            スタティックカラム
メモリアクセスウェイト  なし(040turboの設定)
キャッシュモード        コピーバック
プログラム              gccを使用,-m68040を指定,更に手作業で最適化

●060試作基板
システムクロック        24MHz
メモリモード            スタティックカラム
メモリアクセスウェイト  なし(040turboの設定)
キャッシュモード        コピーバック
プログラム              040turboで使ったものとまったく同じ


○結果

        X68000 PRO              5時間22分
        040turbo                16分3秒
        060試作基板             4分58秒


○補足

 PROの結果は昔記録したもの。
 X68030では最後までテストしていないが1時間以上かかると思われる。
 060モードでは040SYSpatch等の040環境と060spを使用。


                                                        かまだ

060turbo で添付ディスクに入っているプログラムを使うと 2 分 53 秒くらいでした。

使っている画像ファイルは線のパーツを規則的に並べたもの 1 枚だけです。この画像を位置と大きさを固定した SPAN 要素の背景に位置指定付きで貼り付けることで任意の場所に任意の角度のパーツを表示できます。いくつも繋げば直線も曲線も自在に描けるというわけです。長い線を描くには SPAN 要素をいくつも並べなければなりませんが、JavaScript で人工衛星の位置を表示するで使ったように広範囲に渡る線を描くことができます。なお、線ではなくてビットマップを表示したいときはだいぶ前に作ったJavaScript でマンデルブロ集合を描くのように TABLE 要素が使えますが、この方法は大きさに限界があるので人工衛星の軌道の表示には適しません。

東京大学名誉教授で科学雑誌『Newton』を創刊した竹内均さんが亡くなったそうです。毎月読んでいたわけではないのですが、学生時代に『Newton』のフラクタル幾何学の特集を見てフラクタルに興味を持ち始めたことをよく覚えています。掲載されていたフラクタル図形と同じものを MZ-1500 で描いてみたもののどうしても同じ形にならず、結局掲載されていた写真のほうが裏返しだったりもしました。コンピュータグラフィックスの進化した今では珍しくなくなりましたが、『Newton』は創刊当初から比類なき美しさのグラフィックスを用いており、私が頭の体操などで使う画像にちょっと凝りたくなる癖も今思えば『Newton』の影響かも知れません。ご冥福をお祈りいたします。

ニュートン - サイエンス総合情報 (ニュートンプレス)

訃報 (4 月 20 日)

「ニュートン」編集長で地球物理学者の竹内均さん死去 (asahi.com、4 月 20 日)

フラクタル正 8 面体

| コメント(0)

フラクタル正 8 面体の第 6 項です。小さな正 8 面体の数は 46656 個。

このアーカイブについて

このページには、過去に書かれた記事のうちフラクタルカテゴリに属しているものが含まれています。

前のカテゴリは素数です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。

月別 アーカイブ

ウェブページ

Powered by Movable Type 6.3.3