Фрактал Ньютона

2.jpg

Бассейны Ньютона, фракталы Ньютона — разновидность алгебраических фракталов.
Области с фрактальными границами появляются при приближенном нахождении корней нелинейного уравнения алгоритмом Ньютона на комплексной плоскости.

История и дополнительные сведения о фрактале http://ru.wikipedia.org/wiki/Бассейны_Ньютона

Алгоритм
В данном случае рассматривается полином как функция комплексного переменного.
Применим метод Ньютона для нахождения нуля функции комплексного переменного, используя процедуру:
zn+1 = zn - f(zn) / f'(zn)
1. Задается степень полинома, массив его коэффициентов, количество итераций, погрешность.
2. Вычисляется начальное приближение z0 = xmin + ymin * i;
3. Далее в цикле с помощью процедуры описанной выше находятся нули функции,
пока не будет достигнута требуемая погрешность для приближенного корня или заданное количество итераций.

В архиве 2 файла :
complex.c - простая структура для работы с комплексными числами
nf.c - основная программа

#include <stdio.h>
#include <math.h>
#include <conio.h> 
#include <dos.h>
#include <stdlib.h> 
#include <graphics.h>
#include "complex.c"
 
#define EE 1E-80
#define MAX_ROOTS 15
#define BGI_PATH "C:\\Lang\\BCLite\\BGI"
 
typedef struct
{
   Complex root[MAX_ROOTS];
   int nor; // number of root
   int deg; // degree
} Roots;
 
// возвращает номер корня полинома для заполнения области
int process_root(Complex z, Roots *p, float EPS)
{
   int i;
   for(i=0; i < p->nor; i++)
      if (mod( csub( z, p->root[i] ) ) < EPS)
	 return i+1;
   p->root[p->nor]=z;
   p->nor+=1;
   return p->nor;
}
 
// инициализация корней
void Rootinit(Roots *roots)
{
   register int i;
   for(i=0;i<MAX_ROOTS;i++)
      roots->root[i]=comp(0,0);
   roots->nor=0;
   roots->deg=0;
}
 
// возвращает значение полинома в точке z
Complex f(Complex z, int deg, float poly[])
{
   register int i;
   register Complex f;
   f=comp(poly[deg],0);
   for(i=deg-1; i>=0;  i--)
      f=cadd( cmult(f,z), comp(poly[i],0) );
 
   return f;
}
 
// возвращает значение первой производной полинома в точке z
Complex df(Complex z, int deg, float poly[])
{
   register int i;
   register Complex df;
   df=comp(deg*poly[deg],0);
   for(i=deg-1; i>0; i--)
     df=cadd( cmult(df,z), comp(i*poly[i],0) );
 
   return df;
}
 
void Initialize()
{
   int gd, gm, errorcode;
 
   gd = DETECT;
   initgraph(&gd, &gm, BGI_PATH);
 
   errorcode = graphresult();
 
   if (errorcode != grOk)  /* an error occurred */
   {
      printf("Graphics error: %s\n", grapherrormsg(errorcode));
      printf("Press any key to halt:");
      getch();
      exit(1);             /* return with error code */
   }
 
}
int drawFrac(float dx, float per, float poly[], int deg, int NIT, float EPS)
{
    float dy, xmin, ymin, x, y;
    int X,Y;
    int ITEST;
    int i,j,k,kit,col; // счетчики, кол-во итераций, цвет
    Complex z,w,fz,dfz;
    Roots roots;
    FILE *fp;
 
 
    Rootinit(&roots);
    roots.deg = deg;
 
    Initialize();
    per = sqrt(per/100);
    //per *= 0.01;
    // масштаб
    X = getmaxx()*per;
    Y = getmaxy()*per;
    dy= (dx*Y)/X;
    xmin = -dx;
    ymin = dy;
    for(i=0; i<=Y; i++)
    {
       y = ymin - 2*i*dy/Y;
       for(j=0; j<=X; j++)
       {
	 x = xmin + 2*j*dx/X;
	 z=comp(x,y);
	 ITEST=1;
	 // цикл до тех пор, пока не достигнем нужного числа итераций,
	 // или не найдем максимально близкий корень
	 for(k=1; k<=NIT; k++)
	{
	   kit=k;
	   fz = f(z,roots.deg,poly);
	   dfz = df(z,roots.deg,poly);
	   if ( mod(dfz)<=EE ) { ITEST=0; break;}
	   w=z;
	   z=csub( z, cdiv(fz,dfz) ); // zn+1 = zn - f(zn)/f'(z)
	   if( mod( csub(z,w) ) <= EPS) break;
	}
	if (kit == NIT) ITEST=2;
	if (ITEST != 1) putpixel(j,i,BLACK);
	else
	{
	   col = process_root(z,&roots,EPS);
	   putpixel(j,i,col);
	}
      }
      if ( kbhit() != 0) break;
   }
   while(!kbhit()){ }
   restorecrtmode();
   getch();
   return 0;
}
 
 
void main(void)
{
   char answ;
   float poly[MAX_ROOTS+1]={-1,0,0,0,0,1,0};
   float dx=2,per=20;
   float EPS=1E-10;
   int i,deg=5,NIT=40;
 
 
   do
   {
      clrscr();
      printf("INPUT TO DRAW FRACTAL\n");
      printf("1) Polynominal degree (max is %d) : \t%d\n",MAX_ROOTS,deg);
      printf("   Coeficients : \n");
	 for(i=deg;i>=0;i--) printf("\t\tA%d= %f\n",i,poly[i]);
      printf("2) Width of the area on the cartesian surface (dx) : \t%f\n",dx);
      printf("3) Percentage of the area of the screen to be covered: \t%f\n",per);
      printf("4) Maximum Number of Iterations : \t\t\t%d\n",NIT);
      printf("5) Precision : \t\t\t\t\t\t%G \n",EPS);
      printf("6) Draw Fractal\n");
      printf("7) Exit  \n");
      printf("\nWhat is your choice ? ");
      answ=getch();
      clrscr();
      switch (answ)
      {
	case('1'):
	   do
	   {
		 printf("\nEnter the degree of the complex polynominal [<=%d]: ",MAX_ROOTS);
		 scanf("%d",&deg);
	    }while (deg>15 || deg<0);
	    printf("\nEnter the coeficients starting with the greatest power\n");
	    for(i=deg;i>=0;i--)
	    {
		 printf("A%d = ",i);
		 scanf("%f",&poly[i]);
	    }
	    break;
	case('2'):
	   do
	   {
		 printf("\nEnter the width (dx>0) = ");
		 scanf("%f",&dx);
	   }while(dx<=0);
	   break;
	case('3'):
	   do
	   {
		 printf("\nEnter the percentage of the screen to be covered by the fractal(%)\n");
		 printf("(the smallest the fastest ->useful for a preview) : ");
		 scanf("%f",&per);
	   }while (per<=0 || per>100);
	   break;
	case('4'):
	   do
	   {
		 printf("\nEnter the maximum Number of ITerations (NIT>0) : ");
		 scanf("%d",&NIT);
	   }while(NIT<=0);
	   break;
	case('5'):
	    do
	    {
		 printf("\nEnter precision ( >0 ) : ");
		 scanf("%e", &EPS);
	     }while (EPS<0);
	     break;
	case('6'):
	    drawFrac(dx,per,poly,deg,NIT,EPS);
	    break;
	case('7'):
	       exit(0);
     }       			   /*end switch */
  } while(1);                           /*end do-loop*/
 
}

Пример фрактала Ньютона для полинома 5ой степени

Ключевые слова: 
фрактал, бассейны Ньютона
ВложениеРазмер
Nfrac.rar2.58 кб