ニュートン法

ニュートン法でz^3-1=0を解いて収束までの回数と収束先の解が3個存在するうちどれになったかで色分けする。
適当に明るさを割り当てたら明るさが一周して禍々しい…

#include <complex>
#include <cmath>
#include <cstdlib>

#define STB_IMAGE_WRITE_IMPLEMENTATION

#include "stb_image.h"
#include "stb_image_write.h"

using namespace std;

const char* OUT_FILE_NAME = "test.png";

const int N=401;						// 刻み幅の逆数、実数部虚数部それぞれ
									// -1, -1+1/N, -1+2/N, ..., -1+(2N-1)/N, -1+2N/N==1 について調べる
//この数値にあまり意味はなくて、調べてないけどだいたいこれよりいくらか大きくすると画像の出力に失敗する

int main(){
  int width=2*N+1;
  int height=2*N+1;
  int ret;
  char data[2*N+1][2*N+1][3]={};

  complex<double> z,prev,one;
  one.real(1);
  one.imag(0);

  for(int i=0;i<=2*N;++i){
	for(int j=0;j<=2*N;++j){
	  double re=-1.0+(double)i/N, im=-1.0+(double)j/N;
	  z.real(re), z.imag(im);
	  int n;
	  if(0==re&&0==im){
		// 何もしない
	  }else{
		for(n=0;;++n){
		  prev=z;
		  z=(2.0/3.0)*z+(1.0/3.0)*(one/(z*z));
		  if(
			 norm(prev-z) < 1e-12
			 )break;
		}
		data[height-j-1][i][(prev.real()>0.5?0:(prev.imag()>0?1:2))]=255-floor(n*50);
	  }
	}
	// cout << endl;
  }
  ret = stbi_write_png(OUT_FILE_NAME,width,height,3,data,0); //画像の書き出し。この関数は成功すると1を返す
  return !ret;
}

なお画像の書き出しに http://nothings.org/ のライブラリを使った(via http://d.hatena.ne.jp/tueda_wolf/20110702/p2 )。非常に手軽。ただなんか大きい画像を出そうとしたら失敗した(原因不明)。

この場合RGBの順でデータがunsigned charとして並んでおり明るさは0が真っ黒で255が最高っぽい

y軸の方向とかミスっていたので修正

しかしこのへん力学系なのだろうなあ

知り合いは明るさ割り当てるのにatanとか使っていた