質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.50%
C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

Q&A

解決済

1回答

3569閲覧

ニュートン法によるフラクタルが描画できないので、プログラムの何処がダメか指摘して頂きたい。

threeeverytwo

総合スコア49

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

0グッド

0クリップ

投稿2019/01/19 09:59

前提・実現したいこと

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

参考サイトのリンク

http://sky.geocities.jp/bunryu1011/fractalnewton.html

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

episteme

2019/01/19 11:48

「計算できてるけど描画されない」ですか? 「計算できてない」ですか? psetを何度も呼んで線を引く/円を描く はできますか?
threeeverytwo

2019/01/19 12:05

2、3回計算しなおしましたが、計算が出来ていない可能性もあります。
episteme

2019/01/19 12:23

問題の特定はあなたの役目。
threeeverytwo

2019/01/19 12:33

サイトで、x^3+1の解を計算したところ、計算は合っていたので、恐らくコードの部分が問題だと思われます。
threeeverytwo

2019/01/19 12:34

psetを複数回呼び出して、線を引くことは、できました。
episteme

2019/01/19 12:35

なら描画がオカシイんでしょね。
episteme

2019/01/19 12:41

pset(x,y...) してるけど、点(x,y)は描画範囲に収まってますか?
guest

回答1

0

自己解決

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 printf("dx : %f dy : %f\n", dx, dy); 57 58 // ループ処理 59 float x, y, xk, yk; 60 float tmp_x, tmp_y; 61 int cnt_x = 0; 62 int cnt_y = 0; 63 float x2, x3, y2, y3; 64 float a, b, c, d; 65 int flag1, flag2, flag3 = FALSE; 66 67 printf("P1 : %f+%fi P2 : %f+%fi P3 : %f+%fi\n", p1_a, p2_b, p2_a, p2_b, p3_a, p2_b); 68 69 // 窓を開ける 70 screen_int (n+40 , n+40 , RGB (255 ,255 ,255)); 71 //pset(1,1,RGB(255,0,0)); 72 73 74 for ( y = y_min; y <= y_max; y+=dy ) { 75 cnt_x = 0; 76 for ( x = x_min; x <= x_max; x+=dx ) { 77 printf("\nx : %f y : %f\n", x, y); 78 //pset(x, y, RGB(255,0,0)); 79 80 flag1 = FALSE; 81 flag2 = FALSE; 82 flag3 = FALSE; 83 84 tmp_x = x; 85 tmp_y = y; 86 printf("tmp_x : %f tmp_y : %f\n", tmp_x, tmp_y); 87 88 while(1) { 89 x3 = pow(tmp_x, 3.0); 90 x2 = pow(tmp_x, 2.0); 91 y3 = pow(tmp_y, 3.0); 92 y2 = pow(tmp_y, 2.0); 93 94 a = 2 * (x3 - 3 * tmp_x * y2) - 1; 95 b = 2 * (3 * x2 * tmp_y - y3); 96 c = 3 * (x2 - y2); 97 d = 6 * tmp_x * tmp_y; 98 99 if ( pow(c, 2.0) + pow(d, 2.0) == 0 ) { break; } 100 101 xk = (a*c + b*d) / (pow(c, 2.0) + pow(d, 2.0)); 102 yk = (b*c - a*d) / (pow(c, 2.0) + pow(d, 2.0)); 103 104 printf("xk : %f yk : %f\n", xk, yk); 105 //printf("P1 : %f+%fi P2 : %f+%fi P3 : %f+%fi\n", p1_a, p2_b, p2_a, p2_b, p3_a, p2_b); 106 //printf("fabs(xk - p1_a) : %f\n", fabs(xk - p1_a)); 107 //printf("fabs(yk - p1_b) : %f\n", fabs(yk - p1_b)); 108 109 //printf("fabs(xk - p2_a) : %f\n", fabs(xk - p2_a)); 110 //printf("fabs(yk - p2_b) : %f\n", fabs(yk - p2_b)); 111 112 //printf("fabs(xk - p3_a) : %f\n", fabs(xk - p3_a)); 113 //printf("fabs(yk - p3_b) : %f\n\n", fabs(yk - p3_b)); 114 115 if ( fabs(p1_a - xk) <= E ) { 116 printf("ok1\n"); 117 if ( fabsf(p1_b - yk) <= E ) { 118 printf("x : %f y : %f\n", x, y); 119 pset(x+n/2, y+n/2, RGB(255, 0 , 0)); 120 printf("ok11\n"); 121 flag1 = TRUE; 122 } 123 } 124 125 if ( fabsf(p2_a - xk) <= E ) { 126 printf("ok2\n"); 127 if ( fabs(p2_b - yk) <= E ) { 128 printf("x : %f y : %f\n", x, y); 129 pset(x+n/2, y+n/2, RGB(0, 255, 0)); 130 printf("ok22\n"); 131 flag2 = TRUE; 132 } 133 } 134 135 if ( fabsf(p3_a - xk) <= E ) { 136 printf("ok3\n"); 137 if ( fabs(p3_b - yk) <= E ) { 138 printf("x : %f y : %f\n", x, y); 139 pset(x+n/2, y+n/2, RGB(0, 0 , 255)); 140 printf("ok33\n"); 141 flag3 = TRUE; 142 } 143 } 144 145 if ( flag1 || flag2 || flag3 ) { 146 break; 147 } 148 149 //printf("ok4\n"); 150 tmp_x = xk; 151 tmp_y = yk; 152 153 } 154 155 cnt_x++; 156 //printf("cnt_x : %d\n\n", cnt_x); 157 if ( cnt_x == n ) { break; } 158 } 159 cnt_y++; 160 //printf("cnt_y : %d\n\n", cnt_y); 161 if ( cnt_y == n ) { break; } 162 } 163 164 return 0; 165} 166//========================================================================= 167// 関数本体 168//========================================================================= 169 170

投稿2019/01/27 14:13

threeeverytwo

総合スコア49

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

episteme

2019/01/27 21:57

で、ナニがあかんかったですか?
threeeverytwo

2019/01/27 23:55

二重for文内で、x,yにxk,ykを代入していたので、途中でxとyの値が変わっていたのた。つまり、pset()に渡すxとyが毎回同じ値(収束した値)になっていたため。 ただ、x,yに小数点を代入して、何も描画されなかったかは、不明。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.50%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問