Question

Please solve and visualize movement using winbgi.h graphic library ( Visual studio in C program)

•       dx/dt = - y - z

•       dy/dt = x + ay

•       dz/dt = b + z(x - c)

•       Usual parameters: a = b = 0.2, c = 5.7

Visualize z = z(t,x,y).

Using Runga-kutta method and code in "C" language.

This is the information which has been given in the question.

I am looking for Program 'C" code for above question ( visual studio) with mathematical solution.

Answer 1

BEGIN
    REAL c real, c imaginary;
    STRING real and imaginary := IF argc < 3 THEN "-0.8"  ELSE argv( 3 ) FI
                               + " "
                               + IF argc < 4 THEN "0.156" ELSE argv( 4 ) FI
                               + " "
                               ;
    FILE numbers;
    associate( numbers, real and imaginary );
    get( numbers, ( c real, c imaginary ) );
    print( ( fixed( c real, -8, 4 ), fixed( c imaginary, -8, 4 ), newline ) );
    FOR v FROM -100 BY 10 TO 100 DO
        FOR h FROM -280 BY 10 TO 280 DO
            REAL x := h / 200;
            REAL y := v / 100;
            CHAR plot char := "#";
            FOR i TO 50
            WHILE
                REAL z real      = ( x * x ) - ( y * y ) + c real;
                REAL z imaginary = ( x * y * 2 ) + c imaginary;
                IF z real * z real <= 10000
                THEN TRUE
                ELSE
                    plot char := " ";
                    FALSE
                FI
            DO
                x := z real;
                y := z imaginary
            OD;
            print( ( plot char ) )
        OD;
        print( ( newline ) )
    OD
END

Output:

 -0.8000  0.1560



                            # #
                           # #
                        ####   ####
                       ###### ######## #
         ##            ######## ##  # #     ########
            ## #      #########  #      #  ##### # #
        ######## ###    ########     # #    ### #  ## # #
  ####  #####       #     #####     #       #####  ####
# # ##  # ###    # #     ########    ### ########
     # # #####  #      #  #########      # ##
     ########     # #  ## ########            ##
                 # ######## ######
                      ####   ####
                           # #
                          # #

AWK[edit]

Translation of: COBOL

 
# syntax: GAWK -f JULIA_SET.AWK [real imaginary]
BEGIN {
    c_real      = (ARGV[1] != "") ? ARGV[1] : -0.8
    c_imaginary = (ARGV[2] != "") ? ARGV[2] : 0.156
    printf("%s %s\n",c_real,c_imaginary)
    for (v=-100; v<=100; v+=10) {
      for (h=-280; h<=280; h+=10) {
        x = h / 200
        y = v / 100
        plot_char = "#"
        for (i=1; i<=50; i++) {
          z_real = x * x - y * y + c_real
          z_imaginary = x * y * 2 + c_imaginary
          if (z_real ^ 2 > 10000) {
            plot_char = " "
            break
          }
          x = z_real
          y = z_imaginary
        }
        printf("%1s",plot_char)
      }
      printf("\n")
    }
    exit(0)
}
 

Output:

-0.8 0.156



                            # #
                           # #
                        ####   ####
                       ###### ######## #
         ##            ######## ##  # #     ########
            ## #      #########  #      #  ##### # #
        ######## ###    ########     # #    ### #  ## # #
  ####  #####       #     #####     #       #####  ####
# # ##  # ###    # #     ########    ### ########
     # # #####  #      #  #########      # ##
     ########     # #  ## ########            ##
                 # ######## ######
                      ####   ####
                           # #
                          # #

BASIC[edit]

Sinclair ZX81 BASIC[edit]

I don't know exactly how long this takes to run; but I left it for about three and a half hours and when I came back it had already finished. If you can't wait to see the results, I've posted a screenshot here. I also haven't tested it with only 1k of RAM—but I suspect it needs at least 2k.

You can try changing lines 10 and 20 to run the program with different values of the complex constant C+D{\displaystyle i}, or lines 50 and 60 to zoom in.

 10 LET C=-.8
 20 LET D=.156
 30 FOR V=43 TO 0 STEP -1
 40 FOR H=0 TO 63
 50 LET X=(H-32)/21
 60 LET Y=(V-22)/21
 70 FOR A=1 TO 50
 80 LET R=X*X-Y*Y+C
 90 LET I=2*X*Y+D
100 IF R*R>1000 THEN GOTO 150
110 LET X=R
120 LET Y=I
130 NEXT A
140 PLOT H,V
150 NEXT H
160 NEXT V

ZX Spectrum Basic[edit]

Translation of: Sinclair ZX81 BASIC

Higher resolution is obtainable, if you have the time to wait for it.

 10 LET creal=-0.8
 20 LET cimag=0.156
 30 FOR v=-16 TO 16
 40 FOR h=-64 TO 64
 50 LET x=h/40
 60 LET y=v/20
 70 FOR i=1 TO 50
 80 LET zreal=x*x-y*y+creal
 90 LET zimag=x*y*2+cimag
100 IF zreal*zreal>1000 THEN GO TO 150
110 LET x=zreal
120 LET y=zimag
130 NEXT i
140 PLOT h+100,150-v
150 NEXT h
160 NEXT v

Output:

Screenshot here.

C[edit]

Interactive implementation which takes the following 6 parameters as input :

<executable name> <width of graphics window> <height of graphics window> <real part of complex number> <imag part of complex number> <limiting radius> <Number of iterations to be tested>

Prints out usage on incorrect invocation. Requires the WinBGIm library.

 
/*Abhishek Ghosh, 5th October 2017*/
 
#include<graphics.h>
#include<stdlib.h>
#include<math.h>
 
typedef struct{
        double x,y;
}complex;
 
complex add(complex a,complex b){
        complex c;
        c.x = a.x + b.x;
        c.y = a.y + b.y;
        return c;
}
 
complex sqr(complex a){
        complex c;
        c.x = a.x*a.x - a.y*a.y;
        c.y = 2*a.x*a.y;
        return c;
}
 
double mod(complex a){
        return sqrt(a.x*a.x + a.y*a.y);
}
 
complex mapPoint(int width,int height,double radius,int x,int y){
        complex c;
        int l = (width<height)?width:height;
 
        c.x = 2*radius*(x - width/2.0)/l;
        c.y = 2*radius*(y - height/2.0)/l;
 
        return c;
}
 
void juliaSet(int width,int height,complex c,double radius,int n){
        int x,y,i;
        complex z0,z1;
 
        for(x=0;x<=width;x++)
                for(y=0;y<=height;y++){
                        z0 = mapPoint(width,height,radius,x,y);
                        for(i=1;i<=n;i++){
                                z1 = add(sqr(z0),c);
                                if(mod(z1)>radius){
                                        putpixel(x,y,i%15+1);
                                        break;
                                }
                                z0 = z1;
                        }
                        if(i>n)
                                putpixel(x,y,0);
                }
}
 
int main(int argC, char* argV[])
{
        int width, height;
        complex c;
 
        if(argC != 7)
                printf("Usage : %s <width and height of screen, real and imaginary parts of c, limit radius and iterations>");
        else{
                width = atoi(argV[1]);
                height = atoi(argV[2]);
 
                c.x = atof(argV[3]);
                c.y = atof(argV[4]);
 
                initwindow(width,height,"Julia Set");
                juliaSet(width,height,c,atof(argV[5]),atoi(argV[6]));
 
                getch();
        }
 
        return 0;
}
 

C++[edit]

 
#include <windows.h>
#include <string>
#include <complex>
 
const int BMP_SIZE = 600, ITERATIONS = 512;
const long double FCT = 2.85, hFCT = FCT / 2.0;
 
class myBitmap {
public:
    myBitmap() : pen( NULL ), brush( NULL ), clr( 0 ), wid( 1 ) {}
    ~myBitmap() {
        DeleteObject( pen ); DeleteObject( brush );
        DeleteDC( hdc ); DeleteObject( bmp );
    }
    bool create( int w, int h ) {
        BITMAPINFO bi;
        ZeroMemory( &bi, sizeof( bi ) );
        bi.bmiHeader.biSize        = sizeof( bi.bmiHeader );
        bi.bmiHeader.biBitCount    = sizeof( DWORD ) * 8;
        bi.bmiHeader.biCompression = BI_RGB;
        bi.bmiHeader.biPlanes      = 1;
        bi.bmiHeader.biWidth       =  w;
        bi.bmiHeader.biHeight      = -h;
        HDC dc = GetDC( GetConsoleWindow() );
        bmp = CreateDIBSection( dc, &bi, DIB_RGB_COLORS, &pBits, NULL, 0 );
        if( !bmp ) return false;
        hdc = CreateCompatibleDC( dc );
        SelectObject( hdc, bmp );
        ReleaseDC( GetConsoleWindow(), dc );
        width = w; height = h;
        return true;
    }
    void clear( BYTE clr = 0 ) {
        memset( pBits, clr, width * height * sizeof( DWORD ) );
    }
    void setBrushColor( DWORD bClr ) {
        if( brush ) DeleteObject( brush );
        brush = CreateSolidBrush( bClr );
        SelectObject( hdc, brush );
    }
    void setPenColor( DWORD c ) {
        clr = c; createPen();
    }
    void setPenWidth( int w ) {
        wid = w; createPen();
    }
    void saveBitmap( std::string path ) {
        BITMAPFILEHEADER fileheader;
        BITMAPINFO       infoheader;
        BITMAP           bitmap;
        DWORD            wb;
        GetObject( bmp, sizeof( bitmap ), &bitmap );
        DWORD* dwpBits = new DWORD[bitmap.bmWidth * bitmap.bmHeight];
        ZeroMemory( dwpBits, bitmap.bmWidth * bitmap.bmHeight * sizeof( DWORD ) );
        ZeroMemory( &infoheader, sizeof( BITMAPINFO ) );
        ZeroMemory( &fileheader, sizeof( BITMAPFILEHEADER ) );
        infoheader.bmiHeader.biBitCount = sizeof( DWORD ) * 8;
        infoheader.bmiHeader.biCompression = BI_RGB;
        infoheader.bmiHeader.biPlanes = 1;
        infoheader.bmiHeader.biSize = sizeof( infoheader.bmiHeader );
        infoheader.bmiHeader.biHeight = bitmap.bmHeight;
        infoheader.bmiHeader.biWidth = bitmap.bmWidth;
        infoheader.bmiHeader.biSizeImage = bitmap.bmWidth * bitmap.bmHeight * sizeof( DWORD );
        fileheader.bfType    = 0x4D42;
        fileheader.bfOffBits = sizeof( infoheader.bmiHeader ) + sizeof( BITMAPFILEHEADER );
        fileheader.bfSize    = fileheader.bfOffBits + infoheader.bmiHeader.biSizeImage;
        GetDIBits( hdc, bmp, 0, height, ( LPVOID )dwpBits, &infoheader, DIB_RGB_COLORS );
        HANDLE file = CreateFile( path.c_str(), GENERIC_WRITE, 0, NULL, CREATE_ALWAYS, 
            FILE_ATTRIBUTE_NORMAL, NULL );
        WriteFile( file, &fileheader, sizeof( BITMAPFILEHEADER ), &wb, NULL );
        WriteFile( file, &infoheader.bmiHeader, sizeof( infoheader.bmiHeader ), &wb, NULL );
        WriteFile( file, dwpBits, bitmap.bmWidth * bitmap.bmHeight * 4, &wb, NULL );
        CloseHandle( file );
        delete [] dwpBits;
    }
    HDC getDC() const     { return hdc; }
    int getWidth() const  { return width; }
    int getHeight() const { return height; }
    DWORD* bits() const { return ( DWORD* )pBits; }
private:
    void createPen() {
        if( pen ) DeleteObject( pen );
        pen = CreatePen( PS_SOLID, wid, clr );
        SelectObject( hdc, pen );
    }
    HBITMAP bmp; HDC    hdc;
    HPEN    pen; HBRUSH brush;
    void    *pBits; int    width, height, wid;
    DWORD    clr;
};
class julia {
public:
    void draw( std::complex<long double> k ) {
        bmp.create( BMP_SIZE, BMP_SIZE );
        DWORD* bits = bmp.bits();
        int res, pos;
        std::complex<long double> c, factor( FCT / BMP_SIZE, FCT / BMP_SIZE ) ;
 
        for( int y = 0; y < BMP_SIZE; y++ ) {
            pos = y * BMP_SIZE;
 
            c.imag( ( factor.imag() * y ) + -hFCT );
 
            for( int x = 0; x < BMP_SIZE; x++ ) {
                c.real( factor.real() * x + -hFCT );
                res = inSet( c, k );
                if( res ) {
                    int n_res = res % 255;
                    if( res < ( ITERATIONS >> 1 ) ) res = RGB( n_res << 2, n_res << 3, n_res << 4 );
                    else res = RGB( n_res << 4, n_res << 2, n_res << 5 );
                }
                bits[pos++] = res;
            }
        }
        bmp.saveBitmap( "./js.bmp" );
    }
private:
    int inSet( std::complex<long double> z, std::complex<long double> c ) {
        long double dist;//, three = 3.0;
        for( int ec = 0; ec < ITERATIONS; ec++ ) {
            z = z * z; z = z + c;
            dist = ( z.imag() * z.imag() ) + ( z.real() * z.real() );
            if( dist > 3 ) return( ec );
        }
        return 0;
    }
    myBitmap bmp;
};
int main( int argc, char* argv[] ) {
    std::complex<long double> c;
    long double factor = FCT / BMP_SIZE;
    c.imag( ( factor * 184 ) + -1.4 );
    c.real( ( factor * 307 ) + -2.0 );
    julia j; j.draw( c ); return 0;    
}
 

C#[edit]

Translation of: Python

using System.Drawing;
// Note: You have to add the System.Drawing assembly 
//  (right-click "references," Add Reference, Assemblies, Framework,
//    System.Drawing, OK)
using System.Linq;
 
namespace RosettaJuliaSet
{
    class Program
    {
        static void Main(string[] args)
        {
            const int w = 800;
            const int h = 600;
            const int zoom = 1;
            const int maxiter = 255;
            const int moveX = 0;
            const int moveY = 0;
            const double cX = -0.7;
            const double cY = 0.27015;
            double zx, zy, tmp;
            int i;
 
            var colors = (from c in Enumerable.Range(0, 256)
                          select Color.FromArgb((c >> 5) * 36, (c >> 3 & 7) * 36, (c & 3) * 85)).ToArray();
 
            var bitmap = new Bitmap(w, h);
            for (int x = 0; x < w; x++)
            {
                for (int y = 0; y < h; y++)
                {
                    zx = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
                    zy = 1.0 * (y - h / 2) / (0.5 * zoom * h) + moveY;
                    i = maxiter;
                    while (zx * zx + zy * zy < 4 && i > 1)
                    {
                        tmp = zx * zx - zy * zy + cX;
                        zy = 2.0 * zx * zy + cY;
                        zx = tmp;
                        i -= 1;
                    }
                    bitmap.SetPixel(x, y, colors[i]);
                }
            }
            bitmap.Save("julia-set.png");
        }
    }
}
 

C# also makes it relatively easy to do a multi-threaded version, which should run faster than the above:

 
        public struct CalculatedPoint
        {
            public int x;
            public int y;
            public int i;
        }
 
        static void MultiThreaded()
        {
            const int w = 800;
            const int h = 600;
            const int zoom = 1;
            const int maxiter = 255;
            const int moveX = 0;
            const int moveY = 0;
            const double cX = -0.7;
            const double cY = 0.27015;
 
            // Precalculate a pallette of 256 colors
            var colors = (from c in Enumerable.Range(0, 256)
                          select Color.FromArgb((c >> 5) * 36, (c >> 3 & 7) * 36, (c & 3) * 85)).ToArray();
 
            // The "AsParallel" below invokes PLINQ, making evaluation parallel using as many cores as
            // are available.
            var calculatedPoints = Enumerable.Range(0, w * h).AsParallel().Select(xy =>
              {
                  double zx, zy, tmp;
                  int x, y;
                  int i = maxiter;
                  y = xy / w;
                  x = xy % w;
                  zx = 1.5 * (x - w / 2) / (0.5 * zoom * w) + moveX;
                  zy = 1.0 * (y - h / 2) / (0.5 * zoom * h) + moveY;
                  while (zx * zx + zy * zy < 4 && i > 1)
                  {
                      tmp = zx * zx - zy * zy + cX;
                      zy = 2.0 * zx * zy + cY;
                      zx = tmp;
                      i -= 1;
                  }
                  return new CalculatedPoint { x = x, y = y, i = i };
              });
 
            // Bitmap is not multi-threaded, so main thread needs to read in the results as they
            // come in and plot the pixels.
            var bitmap = new Bitmap(w, h);
            foreach (CalculatedPoint cp in calculatedPoints)
                bitmap.SetPixel(cp.x, cp.y, colors[cp.i]);
            bitmap.Save("julia-set-multi.png");
        }