0

C コードでマンデルブロー集合を作成しようとしています。私の出力は、Argandまたは複素平面にプロットされる実部の1列と虚部の1列のデータファイルになります。私はヘッダーで定義されたすべての複雑な数学のものを持っており、実部と虚部に acomplex.hを持つ構造を使用しています。dz のトラフ値をループし、z iter=80 回、またはポイントが定義された半径の外側に来るまで更新しようとしています。dz は虚数であり、本質的に範囲 [(dzrmin + i*dzimin), (dzrmax +i*dzimax)] にあります。z は現在の複素数で、z^2 = dz + z として更新されます。2 つの複素数を合計する関数と、複素数を適切に 2 乗する関数があります。ここに私のコード全体がありますcomplexdouble Rdouble Icsum()csquare()

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include "complex.h"

    int main(void)
    {
        double dzrmin, dzrmax, dzimin, dzimax, dzr,dzi, r0,i0, maxrad; 
        int i,j,k,iter, Ndr,Ndi; 
        complex z, dz; 
        dzrmin = -2.1; 
        dzrmax = 0.6; 
        dzimin = -1.2; 
        dzimax = 1.2; 
        Ndr = 200 ; 
        Ndi = 180; 
        dzr = (dzrmax - dzrmin)/Ndr; 
        dzi = (dzimax - dzimin)/Ndi; 
        r0 = 0.0; 
        i0 = 0.0;
        z.R = r0; 
        z.I = i0; 
        dz.R = dzrmin; 
        dz.I = dzimin; 
        maxrad = 2.0; 
        iter = 80 ; 
        for(i = 0; i < Ndr; i++ )
        {
                    for(j = 0; j < Ndi; j++ )
                    {

                        for(k = 0; k< iter; k++ )
                        {
                    printf("%.6lf %.6lf\n", z.R, z.I ); 
                    z = csum( csquare( z ), dz ) ; 

                    if( cmag(z) > maxrad) k = iter ; 

                }

                dz.I += dzi ; 
            }
            dz.R += dzr; 
            z.R = r0; 
            z.I = i0; 
        }
        return 0; 
    }

そして私のヘッダーファイルcomplex.h

#include <stdlib.h>
#include <math.h>
typedef struct complexnumber{ double R; double I ; } complex ;
double cmag( complex z)                                                                                                                                             
{
    return pow( z.R*z.R + z.I*z.I, 0.5 ) ; 
}

complex csquare( complex z )             //returns square of a complex
{
    complex product ; 
    product.R = z.R*z.R - z.I*z.I ; 
    product.I = 2*z.R*z.I ; 
    return product ; 

}

complex csum( complex z1, complex z2) // sums two complex numbers 
{
    complex sum ; 
    sum.R = z1.R + z2.R ; 
    sum.I = z1.I + z2.I; 
    return sum ; 
}

私はいくつかの実数値を取得しています。それらは非常に急速に大きくなり、半径 2 の外側にあり、多くの実部-nanと虚部が続きます-nan。私が欠けているものに関する提案はありますか?

4

1 に答える 1

1

どのような結果が期待されるかはわかりませんが、このインストルメント化された出力が問題の診断に役立つ可能性があります。

i = 0
i = 0, j = 0
i = 0, j = 0, k = 0: 0.000000 0.000000
i = 0, j = 1
i = 0, j = 1, k = 0: -2.100000 -1.200000
i = 0, j = 2
i = 0, j = 2, k = 0: 0.870000 3.853333
i = 0, j = 3
i = 0, j = 3, k = 0: -16.191278 5.531467
i = 0, j = 4
i = 0, j = 4, k = 0: 229.460353 -180.283027
i = 0, j = 5
i = 0, j = 5, k = 0: 20147.983719 -82736.760384
i = 0, j = 6
i = 0, j = 6, k = 0: -6439430272.999375 -3333957803.416253
i = 0, j = 7
i = 0, j = 7, k = 0: 30350987605860679680.000000 42937577616442236928.000000
i = 0, j = 8
i = 0, j = 8, k = 0: -922453122916892839668491877362832506880.000000 2606395772124638489891541615388582215680.000000
i = 0, j = 9
i = 0, j = 9, k = 0: -5942379156970062175257857153425511563624711669037998408592102022831643293646848.000000 -4808555839107517668369914051798230293111512606121703008400866410836391727464448.000000
i = 0, j = 10
i = 0, j = 10, k = 0: 12189660787377226979177299907109395027773453611835899247334853368367297050334425366457359535808690377953524893091770076537471791519164311649924318416907272192.000000 57148523986878398200502132726020983696472380197135570914789139569106551459309671849901109020634047615447493780748138450365838320365732961897986301575347306496.000000
i = 0, j = 10, k = 1: nan inf
i = 0, j = 10, k = 2: nan nan

内側のループが最初の反復で常に壊れていることに注意してください。問題は、最初から-2.1 - 1.2idzに設定されていることだと思います。しかし、プロットしようとしているマンデルブロ曲線の定義がなければ、何が起こっているのかわかりません。

これは、C99 標準と競合するComplexタイプの代わりに使用して、コードのこの単一ファイル バリアントによって生成されました。complex複雑な関数は静的です。これは、コンパイル時にコンパイラの警告が閉じられるためです。

$ gcc -O3 -g -std=c99 -Wall -Wextra -Wmissing-prototypes -Wstrict-prototypes cx.c -o cx

ソース：

#include <stdio.h>
#include <math.h>

typedef struct Complex{ double R; double I; } Complex;

static double cmag(Complex z)
{
    return pow(z.R*z.R + z.I*z.I, 0.5);
}

static Complex csquare(Complex z)             //returns square of a complex
{
    Complex product;
    product.R = z.R*z.R - z.I*z.I;
    product.I = 2*z.R*z.I;
    return product;
}

static Complex csum(Complex z1, Complex z2) // sums two complex numbers
{
    Complex sum;
    sum.R = z1.R + z2.R;
    sum.I = z1.I + z2.I;
    return sum;
}

int main(void)
{
    double dzrmin, dzrmax, dzimin, dzimax, dzr, dzi, r0, i0, maxrad;
    int i, j, k, iter, Ndr, Ndi;
    Complex z, dz;
    dzrmin = -2.1;
    dzrmax = 0.6;
    dzimin = -1.2;
    dzimax = 1.2;
    Ndr = 200;
    Ndi = 180;
    dzr = (dzrmax - dzrmin)/Ndr;
    dzi = (dzimax - dzimin)/Ndi;
    r0 = 0.0;
    i0 = 0.0;
    z.R = r0;
    z.I = i0;
    dz.R = dzrmin;
    dz.I = dzimin;
    maxrad = 2.0;
    iter = 80;
    for (i = 0; i < Ndr; i++)
    {
        printf("i = %d\n", i);
        for (j = 0; j < Ndi; j++)
        {
            printf("i = %d, j = %d\n", i, j);
            for (k = 0; k < iter; k++)
            {
                printf("i = %d, j = %d, k = %d: ", i, j, k);
                printf("%.6lf %.6lf\n", z.R, z.I);
                z = csum(csquare(z), dz);
                if (cmag(z) > maxrad)
                    break;
            }
            dz.I += dzi;
        }
        dz.R += dzr;
        z.R = r0;
        z.I = i0;
    }
    return 0;
}
于 2013-01-20T01:41:27.387 に答える