global.cpp

V《牛顿法解方程之混沌情况》源代码(C完整应用程序代码)

	double temp5, temp6, temp7, temp8;//, temp9;
	double R1, R2, N7;
	double bx, by;//, xb, yb;
       
	int N, ret;
	int nSum;  //nSum值用来判断迭代是否已经到达根（方程的解），这是一种比较好的方式

	//On Error GoTo aaa:

	x1 = x0;
	y1 = y0;
	N = int(m_dSeData[0][14]);
	N7 = 0.01;
	nSum = 0;
	switch(nN)
	{
	case 1:  //方程 Z^N-1=0
		for( i=1; i<=M; i++)
		{
			P1 = sqrt(x1*x1 + y1*y1);
			SeTa1 = ZArg(x1, y1, 0);
			x2 = ((N-1)*P1*cos(SeTa1) + pow(P1, 1-N)*cos((1-N)*SeTa1))/N;
			y2 = ((N-1)*P1*sin(SeTa1) + pow(P1, 1-N)*sin((1-N)*SeTa1))/N;
			temp0 = pow(fabs(fabs(x1)-fabs(x2)),2) + pow(fabs(fabs(y1)-fabs(y2)),2);

			if(temp0 < N7){
				nSum++;
				if(nSum > 2){
						*dL4 = temp0 / N7;
					if(i > nM){
						break;
					}
				}
			}
			if(i == 1){
				*Hx = fabs(x1 - x2);
				*Hy = fabs(y1 - y2);
				*dL1 = pow((*Hx)*(*Hx)+(*Hy)*(*Hy), 0.5);
			}
			if(i == nM){
				*dL2 = pow((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2), 0.5);
				*dL3 = pow(x2*x2+y2*y2, 0.5);
			}
			x1 = x2;
			y1 = y2;
		}
		ret = (k-1)*M/N+i;
		break;
	case 0: //方程 Z^3-1=0 的特解
	case 2: //方程 Z^3-1=0 的特解
		x1 = -x1;  //其实没有必要，这里做了一下水平翻转
		for(i=1; i<=M; i++)
		{
			x2 = x1;
			y2 = y1;
			m0 = pow((x1*x1 + y1*y1), 2);
			if(pow((x1*x1*x1-x1*y1*y1*3-1),2) + pow((x1*x1*y1*3-y1*y1*y1),2) < N7 )
			{
				*dL4 = fabs(sqrt(x1*x1 + y1*y1) - 1) / N7;
				if(i>nM)
					break;
			}
			x1 = x2 - (x2+y2)/6 + (x2*x2 - y2*y2 - x2*y2*2)/6/m0;
			y1 = y2 + (x2-y2)/6 + (y2*y2 - x2*x2 - x2*y2*2)/6/m0;


			if(i==1){
				*Hx = fabs(x1 - x2);
				*Hy = fabs(y1 - y2);
				*dL1 = pow((*Hx)*(*Hx) + (*Hy)*(*Hy), 0.5);
			}
			if(i==nM){
				*dL2 = pow((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2), 0.5);
				*dL3 = pow(x1*x1 + y1*y1, 0.5);
			}
		}
		if(x1 > 0.9)
			ret = (int)(i + 1 * 0.33 * M);
		else if(y1 > 0.8)
			ret = (int)(i + (0.33 + 1 * 0.33) * M);
		else if(y1 < -0.8)
			ret = (int)(i + (0.66 + 0.33 * 1) * M);
		break;
	case 3:  //方程 Z*(1+Z^A)/(1-Z^A)=R
		bx = x0;
		by = y0;
		A = m_dSeData[0][15]; //: B = SeData(0, 16)
		R1 = m_dSeData[0][17];
		R2 = m_dSeData[0][18];
		for(i=1; i<=M; i++)
		{
			temp1 = bx;
			temp2 = by;
			Zfang(temp1, temp2, A, &temp3, &temp4);
			Zji(temp3, temp4, temp1, temp2, &temp5, &temp6);
			Zji(temp3, temp4, R1, R2, &temp7, &temp8);
			temp5 += temp7 - R1 + temp1;
			temp6 += temp8 - R2 + temp2;
			Zshang(temp3, temp4, temp1, temp2, &temp7, &temp8);
			temp7 += 1 + (A + 1) * temp3;
			temp8 += (A + 1) * temp4;
			Zshang(temp5, temp6, temp7, temp8, &temp3, &temp4);

			bx = temp1 - temp3;
			by = temp2 - temp4;

			temp0 = pow(fabs(fabs(bx)-fabs(temp1)), 2) + 
							pow(fabs(fabs(by)-fabs(temp2)), 2);
			if(temp0 < N7){
				nSum++;
				if(nSum > 2){
					*dL4 = temp0 / N7;
					if(i > nM)
						break;
				}
			}
			if(i==1){
				*Hx = fabs(bx - temp1);
				*Hy = fabs(by - temp2);
				*dL1 = pow(((*Hx)*(*Hx) + (*Hy)*(*Hy)), 0.5);
			}
			if(i==nM){
				*dL2 = pow(((bx-temp1)*(bx-temp1) + (by-temp2)*(by-temp2)), 0.5);
				*dL3 = pow((bx*bx + by*by), 0.5);
			}
		}
		ret = i;
		break;
	case 4:
/*		{
		//    曼德勃罗特集的数学模型非常简单。假设有复数Z和μ，用下式
		//迭代计算：Z=Z2+μ，由于μ的取值不同，经过若干次数迭代以后Z
		//的幅值可能趋向无穷，也可能保持有界，曼德勃罗特集就是那些使
		//Z保持有界的μ的集合，把μ在复平面上的分布作成图像，就像这里
		//演示的那样具有极其复杂的结构。
		int nThre=100;    //nThre代表设置的门限值，当迭代后Z的幅值的平方大于nThre则认为趋于无穷
//		int xCord,yCord;  //xCord和yCord分别代表平面上的点的坐标位置
		int nColor,nTimes;//nColor代表作图时所使用的颜色，nTimes表示迭代的次数
		float reP,rePmin=-1.5,rePmax=0.75; //reP表示由平面上一点所代表的μ值的实部，rePmax、rePmin分别代表μ值实部可取的最大值和最小值
		float imP,imPmin=-1.5,imPmax=1.5; //imP、imPmax和imPmin分别代表μ的虚部及其容许的最大值和最小值
		float reZ,imZ; //reZ和imZ分别代表Z值的实部和虚部
		float deltImP,deltReP;//，deltImP、deltReP表示每一个象素代表的实部和虚部值的大小
		float absZ; //absZ代表Z值的幅值
		float tmpReZ;

		rePmin=(float)G.m_dSeData[0][1];
		imPmin=(float)G.m_dSeData[0][2];
		rePmax=(float)G.m_dSeData[0][3];
		imPmax=(float)G.m_dSeData[0][4];
		//计算屏幕上一个像素表示的实部和虚部值的大小
		deltReP=(rePmax-rePmin)/G.m_dwWidth ;
		deltImP=(imPmax-imPmin)/G.m_dwHeight;
		//计算屏幕上一点所代表的μ值大小
		reP=rePmin+deltReP*x0;
		imP=imPmin+deltImP*y0;
		reZ=0; imZ=0;
		for(nTimes=0;nTimes<M;nTimes++)//160
		{
			tmpReZ=reZ*reZ-imZ*imZ+reP; //计算迭代后的Z
			imZ=2*reZ*imZ+imP;//值及其幅值大小
			reZ=tmpReZ;
			absZ=reZ*reZ+imZ*imZ;
			if (absZ>nThre)
			{
				//nColor=nTimes/10;
				break;
			}
			//如果幅值的平方大于门限值终止迭代
			//nColor=0;
		}
		*Hx = reZ;
		*Hy = imZ;
		*dL1= absZ;
		}*/
		break;
	}
	return ret;
//	Exit Function
//aaa:
//	ret = i;
}

int CGlobal::MDBLT(int x0, int y0)
{
	//    曼德勃罗特集的数学模型非常简单。假设有复数Z和μ，用下式
	//迭代计算：Z=Z2+μ，由于μ的取值不同，经过若干次数迭代以后Z
	//的幅值可能趋向无穷，也可能保持有界，曼德勃罗特集就是那些使
	//Z保持有界的μ的集合，把μ在复平面上的分布作成图像，就像这里
	//演示的那样具有极其复杂的结构。
	//    曼德勃罗特集的数学模型非常简单。假设有复数Z和μ，用下式
	//迭代计算：Z=Z2+μ，由于μ的取值不同，经过若干次数迭代以后Z
	//的幅值可能趋向无穷，也可能保持有界，曼德勃罗特集就是那些使
	//Z保持有界的μ的集合，把μ在复平面上的分布作成图像，就像上
	//面演示的那样具有极其复杂的结构。
	int nThre=10000000;//100;    //nThre代表设置的门限值，当迭代后Z的幅值的平方大于nThre则认为趋于无穷
//		int xCord,yCord;  //xCord和yCord分别代表平面上的点的坐标位置
	int nColor,nTimes;//nColor代表作图时所使用的颜色，nTimes表示迭代的次数
	float reP,rePmin=-1.5,rePmax=0.75; //reP表示由平面上一点所代表的μ值的实部，rePmax、rePmin分别代表μ值实部可取的最大值和最小值
	float imP,imPmin=-1.5,imPmax=1.5; //imP、imPmax和imPmin分别代表μ的虚部及其容许的最大值和最小值
	float reZ,imZ; //reZ和imZ分别代表Z值的实部和虚部
	float deltImP,deltReP;//，deltImP、deltReP表示每一个象素代表的实部和虚部值的大小
	float absZ; //absZ代表Z值的幅值
	float tmpReZ;

	rePmin=(float)G.m_dSeData[0][1];
	imPmin=(float)G.m_dSeData[0][2];
	rePmax=(float)G.m_dSeData[0][3];
	imPmax=(float)G.m_dSeData[0][4];
	//计算屏幕上一个像素表示的实部和虚部值的大小
	deltReP=(rePmax-rePmin)/G.m_dwWidth ;
	deltImP=(imPmax-imPmin)/G.m_dwHeight;
	//计算屏幕上一点所代表的μ值大小
	reP=rePmin+deltReP*x0;
	imP=imPmin+deltImP*y0;
	reZ=0; imZ=0;
	for(nTimes=0;nTimes<256;nTimes++)//160
	{
		tmpReZ=reZ*reZ-imZ*imZ+reP; //计算迭代后的Z
		imZ=2*reZ*imZ+imP;//值及其幅值大小
		reZ=tmpReZ;
		absZ=reZ*reZ+imZ*imZ;
		if (absZ>nThre)
		{
			nColor=nTimes*10;// /10;
			if(nColor>255) nColor%=256;
			break;
		}
		//如果幅值的平方大于门限值终止迭代
		nColor=0;
	}
	return nColor;
	//    上面的程序中，可以通过改变imPmin、imPmax以及rePmin与rePmax的值
	//来对图像的局部细节进行放大，从而可以观察到更加精美复杂的图像，领悟
	//分形图像的神韵。
	//    曼德勃罗特集是人们发现的最早的分形图形之一，也是人们研究最多的
	//分形图形之一，今年的菲尔兹奖（数学界的最高奖）得主麦克马兰就是因为
	//在曼德勃罗特集研究中的成绩而获此殊荣的。这些年关于分形的研究已经渗
	//透到科学领域的各个学科，在计算机领域利用分形的方法来对自然景物进行
	//逼真的模拟是一个很具挑战性的研究方向。
	//      (兰州大学2041信箱 梁昌霖 730000) 
}


//一些复变函数
//============================

//Z1^Z2  （复数的复数次方）
void CGlobal::ZZFang(double x1, double y1, double x2, double y2, double k, double* pX, double* pY)
{
	double T, TT;
	double P, Fai;
//    On Error Resume Next
	P = sqrt(x1 * x1 + y1 * y1);
	if( !P){
		*pX = 0; *pY = 0;
		return;
	}
	Fai = ZArg(x1, y1, k);
	T = pow(P,x2) * exp(-y2 * Fai);
	TT = log(P) * y2 + x2 * Fai;
	*pX = T * cos(TT);
	*pY = T * sin(TT);
}

//Z1*Z2  （复数乘积）
void CGlobal::Zji(double x1, double y1, double x2, double y2, double* pX, double* pY)
{
	*pX = x1*x2 - y1*y2;
	*pY = x1*y2 + y1*x2;
}

//Z1/Z2  （复数商）
void CGlobal::Zshang(double x1, double y1, double x2, double y2, double* pX, double* pY)
{
	double T;
	T = x2*x2 + y2*y2;
	if(!T){
		if( !x1){
			*pX = 1;
			*pY = 1;
		}
		else{
			*pX = (x1>0 ? 1 : (x1<0 ? -1 : 0)) * 1E+50;
			*pY = (y1>0 ? 1 : (y1<0 ? -1 : 0)) * 1E+50;
		}
	}
	else{
		*pX = (x1*x2 + y1*y2) / T;
		*pY = (-x1*y2 + y1*x2) / T;
	}
}

//Z^N  （复数的实数次方）
void CGlobal::Zfang(double x1, double y1, double N, double* pX, double* pY)
{
	double T, TT, AtnYX;
	T = sqrt(x1 * x1 + y1 * y1);
	AtnYX = atan2(y1, x1);
	if(x1 < 0)
		TT = PI + AtnYX;
	else if(y1 > 0)
		TT = AtnYX;
	else //if y1<=0 then
		TT = PI * 2 + AtnYX;//2#

	T = pow(T,N);
	TT = TT * N;
	*pX = T * cos(TT);
	*pY = T * sin(TT);
}

//Arg(Z)  （复数的辐角）
double CGlobal::ZArg(double x, double y, double k)
{
	double ret=0;
	if(!x){
		if(!y)
			ret = 0;
		else if(y > 0)
			ret = PI / 2;//2#
		else //if y<0 then
			ret = PI * 1.5;
	}
	else if(x>0){
		ret = atan2(y, x);
		if(ret < 0)
			ret = PI * 2 + ret;
	}
	else //if x<0 then
		ret = atan2(y, x) + PI;

	ret = ret + PI * 2 * k;

	return ret;
}

//为了实现 Mandelbrot(曼德勃罗特集) 特效定义的函数 （其实就是Mandelbrot函数迭代）
void CGlobal::fz2(double x0, double y0, double xx, double yy, double* pX, double* pY, int N)
{//Z=Z2+μ
	double t1, t2, t3, t4;
	int i;

	t1 = x0;
	t2 = y0;
	for(i=1; i<=N; i++){