1500字范文,内容丰富有趣,写作好帮手!
1500字范文 > 数学分形之复牛顿迭代法

数学分形之复牛顿迭代法

时间:2019-02-20 23:29:35

相关推荐

数学分形之复牛顿迭代法

复牛顿迭代法

绘图库:Easy Graphics Engine (EGE)

编程语言:c++

z=z^3-1

z=z^4-1

z=sin(z)

z=z^6-1

代码:

#include <stdio.h>#include <stdlib.h>#include <math.h>#include <string.h>#include <windows.h>#include <graphics.h>#include <conio.h>#include <graphics.h>#include <complex.h>#define pi 3.1415926#define H 600.0#define W 600.0const int NUM_ITER = 1000;// 迭代次数const int N = 500;#ifndef SWAP#define SWAP(a, b, t) {t = a; a = b; b = t;}#endif // ! SWAPvoid color(short x)//设置字体颜色的自定义函数{if((x>=0)&&(x<=15))SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),x);elseSetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),7);}//0=黑色 1=蓝色 2=绿色 3=湖蓝色 4=红色 5=紫色//6=黄色 7=白色 8=灰色 9=淡蓝色 10=淡绿色 11=淡浅绿色//12=淡红色 13=淡紫色 14=淡黄色 15=亮白色struct Complex{double re, im;Complex() :re(0.0), im(0.0) {}Complex(double real, double imag) : re(real), im(imag) {}//重载运算符Complex operator * (Complex c) {return Complex(re * c.re - im * c.im, im * c.re + re * c.im); }//Complex operator / (Complex c) { if(c.re==0&&c.im==0){return Complex(0,0);}else {return Complex((re*c.re+im*c.im)/(c.re*c.re+c.im*c.im),(im*c.re-re*c.im)/(c.re*c.re+c.im*c.im))}};Complex operator + (Complex c) {return Complex(re + c.re, im + c.im); }Complex operator - (Complex c) {return Complex(re - c.re, im - c.im); }//Complex operator ^ (double n) { return Complex(pow(sqrt(re*re+im*im),n)*cos(n*atan(im/re)),pow(sqrt(re*re+im*im),n)*sin(n*atan(im/re))); }};Complex operator / (Complex a, Complex b){Complex c;if(b.re==0&&b.im==0){c.re = 0;c.im = 0;}else{c.re=(a.re*b.re+a.im*b.im)/(b.re*b.re+b.im*b.im);c.im=(a.im*b.re-a.re*b.im)/(b.re*b.re+b.im*b.im);}return c;}Complex operator ^ (Complex a,double n){Complex c;double r=sqrt(a.re*a.re+a.im*a.im);if(a.re!=0){c.re=pow(r,n)*cos(n*atan(a.im/a.re));c.im=pow(r,n)*sin(n*atan(a.im/a.re));}else if(a.im*a.re>0){c.re=pow(r,n)*cos(n*(pi/2));c.im=pow(r,n)*sin(n*(pi/2));}else{c.re=pow(r,n)*cos(n*(-pi/2));c.im=pow(r,n)*sin(n*(-pi/2));}return c;}// 初始化颜色#define MAX_COLOR_NUM 64int Color[MAX_COLOR_NUM];void initColdeTable(){for (int i = 0; i < MAX_COLOR_NUM / 2; i++){Color[i] = HSLtoRGB(360, 1.0, i * 2.0 / MAX_COLOR_NUM);Color[MAX_COLOR_NUM - 1 - i] = HSLtoRGB(0, 1.0, i * 2.0 / MAX_COLOR_NUM);}}//putpixel(x,y, (i == NUM_ITER) ? 0 : Color[i % MAX_COLOR_NUM]);// 绘制曼德布洛特集 (Mandelbrot Set)void draw(double fromx, double fromy, double tox, double toy){Complex z,zz,e,t;t.re=7.0;t.im=0;e.re=1.0;e.im=0;double r;for (int x = 0; x < W; x++){z.re = fromx + (tox - fromx)*(x/W);for (int y = 0; y < H; y++){z.re = fromx + (tox - fromx)*(x/W);z.im=fromy+(toy-fromy)*(y/H);int i;r=1;for (i=0;i<100; i++){zz.re=z.re;zz.im=z.im;z=z-(z*z*z*z*z*z*z-e)/(t*z*z*z*z*z*z);r=sqrt((z.re-zz.re)*(z.re-zz.re)+(z.im-zz.im)*(z.im-zz.im));if(r<0.0000001) break;}//putpixel(x,-y+600,Color[i%MAX_COLOR_NUM]);putpixel(x,y, (i == NUM_ITER) ? 0 : Color[i % MAX_COLOR_NUM]);//putpixel(x, y,RGB(20*log(i),20*log(i),20*log(i)));}}}// 绘制茱莉亚集void draw1(double fromx, double fromy, double tox, double toy,double xx,double yy){Complex z,c;c.re = xx, c.im = yy;int k;for (int x = 0; x < W; x++){for (int y = 0; y < H; y++){z.re = fromx + (tox - fromx) * (x / W);z.im = fromy + (toy - fromy) * (y / H);int i = 0;for ( ; i < NUM_ITER; i++){if (z.re * z.re + z.im * z.im > 4.0)break;z = z * z + c;}putpixel(x, y, (i == NUM_ITER) ? 0 : Color[i % MAX_COLOR_NUM]);}}setbkmode(TRANSPARENT);setcolor(WHITE);setfont(20, 10, "黑体");outtextxy(250, 30, "朱利亚集合");setfont(20, 10, "黑体");outtextxy(250, 50, "Julia set");setfont(20,10, "黑体");outtextxy(400, 550, "c.real :");outtextxy(400, 570, "c.image:");char s[100],f[100];sprintf(s, "%lf", xx);sprintf(f, "%lf", yy);outtextxy(500, 550,s);outtextxy(500, 570,f);}int main(){double x,y;color(11);//printf(" 茱莉亚集\n\n");ccc:initgraph(W,H, INIT_RENDERMANUAL);/* color(11);printf("请输入复数c的实部a,虚部b\n");color(12);scanf("%lf %lf",&x,&y);color(10);printf("按Esc键重置程序");printf("\n\n\n");*///初始化颜色表initColdeTable();// 初始化 Mandelbrot Set(曼德布洛特集)坐标系范围const double INIT_FROM_X = -2, INIT_TO_X = 2;const double INIT_FROM_Y = 2, INIT_TO_Y= -2;double fromx = INIT_FROM_X, tox = INIT_TO_X;double fromy = INIT_FROM_Y, toy = INIT_TO_Y;draw(fromx, fromy, tox, toy);printf("绘制完毕\n\n");bool isPress = false;//鼠标左键按下标志位bool redraw = true;int areaLeft = 0, areaTop = 0, areaRight = 0, areaBottom = 0;// 定义选区while (1){mouse_msg msg = getmouse();// 鼠标中键按下时重置图形if (msg.is_mid() && msg.is_down()){fromx = INIT_FROM_X;tox = INIT_TO_X;fromy = INIT_FROM_Y;toy = INIT_TO_Y;redraw = true;}//鼠标左键点击,选取范围else if (msg.is_left()) {if (msg.is_down()) {isPress = true;setcolor(WHITE);setwritemode(R2_XORPEN);areaLeft = areaRight = msg.x;areaTop = areaBottom = msg.y;}else {// 鼠标左键松开时确定选区isPress = false;redraw = true;//消除选框rectangle(areaLeft, areaTop, areaRight, areaBottom);setwritemode(R2_COPYPEN);areaRight = msg.x;areaBottom = msg.y;if (areaLeft != areaRight && areaTop != areaBottom) {// 修正选区为 4:4int temp;if (areaLeft > areaRight)SWAP(areaLeft, areaRight, temp);if (areaTop > areaBottom)SWAP(areaTop, areaBottom, temp);if ((areaRight - areaLeft) * 0.75 < (areaBottom - areaTop)){areaBottom += (3 - (areaBottom - areaTop) % 3);areaLeft -= (areaBottom - areaTop) / 4 * 4 / 2 - (areaRight - areaLeft) / 2;areaRight = areaLeft + (areaBottom - areaTop) / 4 * 4;}else{areaRight += (4 - (areaRight - areaLeft) % 4);areaTop -= (areaRight - areaLeft) * 4 / 4 / 2 - (areaBottom - areaTop) / 2;areaBottom = areaTop + (areaRight - areaLeft) * 4 / 4;}// 更新坐标系double from = fromx, to = tox;fromx = from + (to - from) * areaLeft / W;tox = from + (tox - from) * areaRight / W;from = fromy;to = toy;fromy = from + (to - from) * areaTop / H;toy = from + (to - from) * areaBottom /H;}}}else if (msg.is_move() && isPress){//消除选框rectangle(areaLeft, areaTop, areaRight, areaBottom);areaRight = msg.x;areaBottom = msg.y;//绘制选框rectangle(areaLeft, areaTop, areaRight, areaBottom);}//重绘if (redraw){redraw = false;draw(fromx, fromy, tox, toy);}if (kbhit()){if (getch() == 27) goto ccc;}}getch();closegraph();return 0;}

本内容不代表本网观点和政治立场,如有侵犯你的权益请联系我们处理。
网友评论
网友评论仅供其表达个人看法,并不表明网站立场。