前提・実現したいこと
z^3 + 1 = 0 をニュートン法で解く場合の初期値設定問題を視覚的に理解するた
めに,引力圏図(basin of attraction)をC 言語を用いてフラクタルを描く。
ニュートン法の公式
a1 = a0 - f(a0)/f'(a0) … (1)
a1の値を(1)のa0に代入すると、a2の値が得られるはず。
そして、次のa2の値を(1)に代入すると、a3の値が得られる。
これを繰り返すと、ある値に収束するはず。
z^3 + 1 = 0 より
z = -1 , (1+√3i)/2, (1-√3i)/2 なので
この3つの値に収束するはずである。
z = x + bi を代入して、実数と虚数に分けてニュートン法を用い、精度0.01で収束したら、点を描画する
という、プログラムなのだが、全く描画されない... (T T)。
試したこと
配布されている資料やサイトを閲覧したり、計算を紙で解いてみたり、printして、数値を見ながらデバッグしたり、考えつく限りの方法を試した。
該当のソースコード
c
1//========================================================================= 2// ライブラリ 3//========================================================================= 4#include <stdio.h> 5#include <math.h> 6#pragma comment (lib , "user32.lib") 7#pragma comment (lib , "gdi32.lib") 8#include <windows.h> 9#include "screen_int.f" 10#include "pset.f" 11//========================================================================= 12// 定義 13//========================================================================= 14#define E 0.01 15#define TRUE 1 16#define FALSE 0 17//========================================================================= 18// 関数定義 19//========================================================================= 20 21//========================================================================= 22// 本体関数 23//========================================================================= 24int main(void) { 25 26 // 初期宣言 27 float x_min, x_max; 28 float y_min, y_max; 29 int n; 30 float p1_a = -1; 31 float p1_b = 0; 32 float p2_a = 0.5; 33 float p2_b = -1*sqrtf(3)/2; 34 float p3_a = 0.5; 35 float p3_b = sqrtf(3)/2; 36 37 // 入力処理 38 printf("x_min:"); 39 scanf("%f", &x_min); 40 41 printf("x_max:"); 42 scanf("%f", &x_max); 43 44 printf("y_min:"); 45 scanf("%f", &y_min); 46 47 printf("y_max:"); 48 scanf("%f", &y_max); 49 50 printf("n :"); 51 scanf("%d", &n); 52 53 float dx = (x_max - x_min) / n; 54 float dy = (y_max - y_min) / n; 55 56 // ループ処理 57 float x, y, xk, yk; 58 float tmp_x, tmp_y; 59 int cnt_x = 0; 60 int cnt_y = 0; 61 float x2, x3, y2, y3; 62 float a, b, c, d; 63 int flag1, flag2, flag3 = FALSE; 64 65 printf("P1 : %f+%fi P2 : %f+%fi P3 : %f+%fi\n", p1_a, p2_b, p2_a, p2_b, p3_a, p2_b); 66 67 // 窓を開ける 68 screen_int (n+40 , n+40 , RGB (255 ,255 ,255)); 69 70 for ( y = y_min; y <= y_max; y+=dy ) { 71 cnt_x = 0; 72 for ( x = x_min; x <= x_max; x+=dx ) { 73 printf("x : %f y : %f\n", x, y); 74 flag1 = FALSE; 75 flag2 = FALSE; 76 flag3 = FALSE; 77 tmp_x = x; 78 tmp_y = y; 79 80 while(1) { 81 x3 = pow(x, 3.0); 82 x2 = pow(x, 2.0); 83 y3 = pow(y, 3.0); 84 y2 = pow(y, 2.0); 85 86 a = 2 * (x3 - 3 * x * y2) - 1; 87 b = 2 * (3 * x2 * y - y3); 88 c = 3 * (x2 - y2); 89 d = 6 * x * y; 90 91 xk = (a*c + b*d) / (pow(c, 2.0) + pow(d, 2.0)); 92 yk = (b*c - a*d) / (pow(c, 2.0) + pow(d, 2.0)); 93 94 printf("xk : %f yk : %f\n", xk, yk); 95 printf("P1 : %f+%fi P2 : %f+%fi P3 : %f+%fi\n", p1_a, p2_b, p2_a, p2_b, p3_a, p2_b); 96 printf("fabs(xk - p1_a) : %f\n", fabs(xk - p1_a)); 97 printf("fabs(yk - p1_b) : %f\n", fabs(yk - p1_b)); 98 99 printf("fabs(xk - p2_a) : %f\n", fabs(xk - p2_a)); 100 printf("fabs(yk - p2_b) : %f\n", fabs(yk - p2_b)); 101 102 printf("fabs(xk - p3_a) : %f\n", fabs(xk - p3_a)); 103 printf("fabs(yk - p3_b) : %f\n\n", fabs(yk - p3_b)); 104 105 if ( fabs(xk - p1_a) <= E ) { 106 printf("ok1\n"); 107 if ( fabsf(yk - p1_b) <= E ) { 108 pset(x, y, RGB(255, 0 , 0)); 109 printf("ok11\n"); 110 flag1 = TRUE; 111 } 112 } 113 114 if ( fabsf(xk - p2_a) <= E ) { 115 printf("ok2\n"); 116 if ( fabs(yk - p2_b) <= E ) { 117 pset(x, y, RGB(0, 255 , 0)); 118 printf("ok22\n"); 119 flag2 = TRUE; 120 } 121 } 122 123 if ( fabsf(xk - p3_a) <= E ) { 124 printf("ok3\n"); 125 if ( fabs(yk - p3_b) <= E ) { 126 pset(x, y, RGB(0, 0 , 255)); 127 printf("ok33\n"); 128 flag3 = TRUE; 129 } 130 } 131 132 if ( flag1 || flag2 || flag3 ) { 133 break; 134 } 135 136 printf("ok4\n"); 137 x = xk; 138 y = yk; 139 140 } 141 x = tmp_x; 142 y = tmp_y; 143 cnt_x++; 144 printf("cnt_x : %d\n\n", cnt_x); 145 if ( cnt_x == n ) { break; } 146 } 147 cnt_y++; 148 printf("cnt_y : %d\n\n", cnt_y); 149 if ( cnt_y == n ) { break; } 150 } 151 152 return 0; 153} 154//========================================================================= 155// 関数本体 156//========================================================================= 157 158
参考サイトのリンク
回答1件
あなたの回答
tips
プレビュー