0
点赞
收藏
分享

微信扫一扫

数学分形之复牛顿迭代法

Silence潇湘夜雨 2022-01-17 阅读 72
c++几何学数学

复牛顿迭代法

绘图库：Easy Graphics Engine （EGE）
编程语言：c++
z=z^3-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.0


const int NUM_ITER = 1000;		// 迭代次数
const int N = 500;
#ifndef  SWAP
#define SWAP(a, b, t) {t = a; a = b; b = t;}
#endif // ! SWAP
void color(short x)//设置字体颜色的自定义函数
{
    if((x>=0)&&(x<=15))
        SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE),x);
    else
        SetConsoleTextAttribute(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 64
int 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:4
					int 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;
}
举报

相关推荐

不动点迭代法和牛顿迭代法

python迭代牛顿法迭代法Python编程语言

牛顿迭代法leetcode

leetcode算法职场和发展

牛顿迭代法求开方

牛顿迭代法计算机编程数值分析OpenStack云计算

c++:牛顿迭代法

c++

梯度下降法，牛顿法，高斯-牛顿迭代法

线性代数c++

Logistic回归与牛顿迭代法

i++Data牛顿迭代法HarmonyOS后端开发

牛顿迭代法求平方根

算法动态规划几何学

YTU 2405: C语言习题 牛顿迭代法求根

主函数ios牛顿迭代法虚拟化云计算

java迭代法

迭代迭代法Java后端开发

迭代法求解ICP

算法公式推导github.netPython后端开发
0 条评论