Trouble with creating Mandelbrot Set plot in C

Question

I am trying to create a Mandelbrot set in C code; my output will be a data file, one column of the real parts and one column of the imaginary parts to be plotted in the Argand or complex plane. I have all the complex math stuff defined in a header complex.hand am using a structure complex with a double R for the real part and double I for the imaginary. I am trying to loop trough values for dz and update z iter=80 number of times, or until the point lies outside a defined radius. dz is imaginary, essentially on a range [(dzrmin + i*dzimin), (dzrmax +i*dzimax)]. z is the current complex number, updated as z^2 = dz + z. I have a function csum() to sum two complex numbers and a function csquare() to properly square a complex number. Here is my whole code

    #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; 
    }

and my header file 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 ; 
}

I am getting a few real values, they get very big really quickly and lie outside a radius of 2, followed by a lot of real parts -nan and imaginary parts -nan. Any suggestions as to what I am missing?

Answer 1

I'm not sure what results you're expecting, but maybe this instrumented output will help you diagnose the problem:

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

Note how the inner loop is breaking on the first iteration, all the time. I think the trouble is that dz is set to -2.1 - 1.2i to start with. But without a definition of the Mandelbrot curve you're seeking to plot, I can't tell what should be going on.

This was produced by this single-file variant of your code, using Complex in place of the type complex which conflicts with the C99 standard. The complex functions are static because that shuts up the compiler warnings when I compile with:

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

Source:

#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;
}