equation.cpp

求解方程根的函数和一些微分方程迭代的函数

#include "equation.h"
#include "math.h"
#include "iostream.h"
#include "matrix.h"

int equat_newton(double *a,double *xr,double *xi,int n)
{
	//求实系数代数方程全部根的牛顿-下山法
	//其中a(n+1),存放多项式方程的实系数(从后面开始放)
    //n整形变量，多项式方程的次数
    //xr返回n个根的实部，xr返回n个根的虚部
	int m,i,jt,k,is,it;
	double t,x,y,x1,y1,dx,dy,p,q,w,dd,dc,c;
	double g,u,v,pq,g1,u1,v1;
	void g60(double *,double *,double *,double *,double*,double *,double *,
		double *,double *,int *,int *);
	void g65(double *,double *,double *,double *,double *,double *,double *,
		double *,double *,int *,int *,int *);
	void g90(double xr[],double xi[],double *,double *,double *,double *,
		double *,double *,int *);
	m=n;
	while((m>0)&&(fabs(a[m])+1.0==1.0))
		m=m-1;  //要求全部系数都不等于零
	if(m<=0)
	{
		cout<<"fail!\n";
		return (-1);
	}
	for(i=0;i<=m;i++)
		*(a+i)=*(a+i)/(0+*(a+m));
	for(i=0;i<=m/2;i++)
	{
		w=*(a+i);
		*(a+i)=*(a+m-i);
		*(a+m-i)=w;
	}
	k=m;
	is=0;
	w=1.0;
	jt=1;
	while(jt==1)
	{
		pq=fabs(*(a+k));
		while(pq<1.0e-12)
		{
			//该while语句当系数小于一定值，当作零
			xr[k-1]=0.0;
			xi[k-1]=0.0;
			k=k-1;
			if(k==1)
			{
				xr[0]=-*(a+1)*w/(*a);
				xi[0]=0.0;
			}
			pq=fabs(*(a+k));
		}
		q=log(pq);
		q=q/(1.0*k);
		q=exp(q);
		p=q;
		w=w*p;
		for(i=1;i<=k;i++)
		{
			*(a+i)=*(a+i)/q;
			q=q*p;
		}
		x=0.0001;
		x1=x;
		y=0.2;
		y1=y;
		dx=1.0;
		g=1.0e+37;
		biaoqian:
		u=*a;
		v=0.0;
		for(i=1;i<=k;i++)
		{
			p=u*x1;
			q=v*y1;
			pq=(u+v)*(x1+y1);
			u=p-q+*(a+i);
			v=pq-p-q;
		}
		g1=u*u+v*v;
		if(g1>=g)
		{
			if(is!=0)
			{
				it=1;
				g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
				if(it==0)  goto biaoqian;
			}
			else
			{
				g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
				if(t>=1.0e-3) goto biaoqian;
				if(g>1.0e-18)
				{
					it=0;
					g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
					if(it==0)  goto biaoqian;
				}
			}
			g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
		}
		else
		{
			g=g1;
			x=x1;
			y=y1;
			is=0;
			if(g<=1.0e-22)
				g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
			else
			{
				u1=k*(*a);
				v1=0.0;
				for(i=2;i<=k;i++)
				{
					p=u1*x;
					q=v1*y;
					pq=(u1+v1)*(x+y);
					u1=p-q+(k-i+1)**(a+i-1);
					v1=pq-p-q;
				}
				p=u1*u1+v1*v1;
				if(p<=1.0e-20)
				{
					it=0;
					g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
					if(it==0)  goto biaoqian;
					g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
				}
				else
				{
					dx=(u*u1+v*v1)/p;
					dy=(u1*v-v1*u)/p;
					t=1.0+4.0/k;
					g60(&t,&x,&y,&x1,&y1,&dx,&dy,&p,&q,&k,&it);
					if(t>=1.0e-03) goto biaoqian;
					if(g>1.0e-18)
					{
						it=0;
						g65(&x,&y,&x1,&y1,&dx,&dy,&dd,&dc,&c,&k,&is,&it);
						if(it==0) goto biaoqian;
					}
					g90(xr,xi,a,&x,&y,&p,&q,&w,&k);
				}
			}
		}
		if(k==1)
			jt=0;
		else
			jt=1;
	}
		return 1;
}

void g60(double *t,double *x,double *y,double *x1,double *y1,double *dx,
		 double *dy,double *p,double *q,int *k,int *it)
{
	//牛顿-下山法求解方程的子程序
	*it=1;
	while(*it==1)
	{
		*t=*t/1.67;
		*it=0;
		*x1=*x-(*t)*(*dx);
		*y1=*y-(*t)*(*dy);
		if(*k>=50)
		{
			*p=sqrt((*(x1)*(*x1)+*(y1)*(*y1)));
			*q=exp(85.0/(*k));
			if(*p>=*q)
				*it=1;
		}
	}
	return ;
}
void g90(double xr[],double xi[],double *a,double *x,double *y,double *p,
		 double *q,double *w,int *k)
{
	//牛顿-下山法求解方程的子程序
	int i;
	if(fabs(*y)<=1.0e-06)
	{
		*p=-(*x);
		*y=0.0;
		*q=0.0;
	}
	else
	{
		*p=-2.0*(*x);
		*q=(*x)*(*x)+(*y)*(*y);
		xr[*k-1]=(*x)*(*w);
		xi[*k-1]=-(*y)*(*w);
		*k=*k-1;
	}
	for(i=1;i<=*k;i++)
	{
		*(a+i)=*(a+i)-*(a+i-1)*(*p);
		*(a+i+1)=*(a+i+1)-*(a+i-1)*(*q);
	}
	xr[*k-1]=(*x)*(*w);
	xi[*k-1]=(*y)*(*w);
	*k=*k-1;
	if(*k==1)
	{
		xr[0]=-*(a+1)*(*w)/(*a);
		xi[0]=0.0;
	}
	return ;
}

void g65(double *x,double *y,double *x1,double *y1,double *dx,double *dy,
		 double *dd,double *dc,double *c,int *k,int *is,int *it)
{
	//牛顿-下山法求解方程的子程序
	if(*it==0)
	{
		*is=1;
		*dd=sqrt((*dx)*(*dx)+(*dy)*(*dy));
		if(*dd>1.0)  
			*dd=1.0;
		*dc=6.28/(4.5*(*k));
		*c=0.0;
	}

	while(1==1)
	{
		*c=*c+*dc;
		*dx=(*dd)*cos(*c);
		*dy=(*dd)*sin(*c);
		*x1=*x+*dx;
		*y1=*y+*dy;
		if(*c<=6.29)
		{
			*it=0;
			return ;
		}
		*dd=*dd/1.67;
		if(*dd<=1.0e-07)
		{
			*it=1;
			return ;
		}
		*c=0.0;
	}
}

void equat_rungekutta(double t,double h,double *y,double *z,int n,int k)
{
	//全区间积分的定步长龙格-库塔法
    //t为积分的起始点，h为定步长，y存放初始值
    //z返回k个积分点上的n个未知函数值，n微分函数个数，k积分的步数（含积分点）
	extern void rungekuttaf(double t,double *y,double *d,int n);
	int i,j,l;
	double a[4],tt,b[20],d[20];//b存放y值，d存放微分函数得到的新y值？
	a[0]=h/2;
	a[1]=h/2;
	a[2]=h;
	a[3]=h;
	for(i=0;i<n;i++)
		z[i*k]=y[i];//得到第一组的y值
	for(l=1;l<k;l++)//k组数据
	{
		rungekuttaf(t,y,(double *)d,n);//需要自己编辑微分方程
		for(i=0;i<n;i++)
			b[i]=y[i];
		for(j=0;j<3;j++)
		{
			for(i=0;i<n;i++)
			{
				y[i]=y[i]+a[j]*d[i];
				b[i]=b[i]+a[j+1]*d[i]/3.0;
			}
			tt=t+a[j];
			rungekuttaf(tt,y,(double *)d,n);
		}
		for(i=0;i<n;i++)
			y[i]=b[i]+h*d[i]/6.0;
		for(i=0;i<n;i++)
			z[i*k+l]=y[i];
		t=t+h;
	}
	return ;
}

void equat_rungekutta(double t,double h,double *y,int n,double eps)
{
	//积分一步的变步长龙格-库塔法(与前面函数重载)
    //t为积分的起始点，h为定步长，y存放初始值,返回时存放结果
    //n微分函数个数，eps积分精度要求。 
    //每次调用该函数的时候同样需要自己写一个微分方程的函数
	extern void rungekuttaf(double t,double *y,double *d,int n);
	int m,i,j,k;
	double hh,p,dt,x,tt,q,a[4];
	double g[20],b[20],c[20],d[20];
	hh=h;
	m=1;//用做步长变为一半
	p=1.0+eps;
	x=t;
	for(i=0;i<n;i++)
	{
		for(i=0;i<n;i++)
			c[i]=y[i];//保存初值
		while(p>=eps)
		{
			a[0]=hh/2;
		    a[1]=hh/2;
		    a[2]=hh;
		    a[3]=hh;
			for(i=0;i<n;i++)
			{
				g[i]=y[i];//保存计算得到的一步结果
				y[i]=c[i];//重新导入初始值
			}
			dt=h/m;
			t=x;
			for(j=0;j<m;j++)  //当m变大的时候，就变成了全区间积分了
			{
				rungekuttaf(t,y,(double *)d,n);
				for(i=0;i<n;i++)
					b[i]=y[i];//与前面算法一致
				for(k=0;k<3;k++)
				{
					for(i=0;i<n;i++)
					{
						y[i]=y[i]+a[k]*d[i];
						b[i]=b[i]+a[k+1]*d[i]/3.0;
					}
					tt=t+a[k];
					rungekuttaf(tt,y,(double *)d,n);
				}
				for(i=0;i<n;i++)
					y[i]=b[i]+hh*d[i]/6.0;
				t=t+dt;
			}
			p=0.0;
			for(i=0;i<n;i++)  //检查精度要求
			{
				q=fabs(y[i]-g[i]);
				if(q>p) 
					p=q;
			}
			hh=hh/2.0;
			m=m+m;
		}
	}
	return ;
}
				
/*void equat_gill(double t,double h,double eps,double *y,double *q,int n)
{
	extern void gillf(double t,double *y,double *d,int n);
	int i,j,k,m,ii;
	double x,p,hh,r,s,t0,dt,qq;
	double d[20],u[20],v[20],g[20];
	//求y[i]的一些参数
	double a[4]={0.5,0.29289321881,1.7071067812,0.166666667};
	double b[4]={2.0,1.0,1.0,2.0};//y[i]中q的系数
	double c[4];//q[i]中k的系数(1)
	double e[4]={0.5,0.5,1.0,1.0};//q[i]中k的系数(2)
	for(i=0;i<3;i++)
		c[i]=a[i];
	c[3]=0.5;
	x=t;
	p=1.0+eps;
	hh=h;
	m=1;
	for(j=0;j<n;j++)
		u[j]=*(y+j);  //保存初始值
	while(p>=eps)
	{   for(j=0;j<n;j++)
		{
			v[j]=*(y+j);
		    *(y+j)=u[j];  //保存初始值
		    g[j]=q[j];  
		}
	     dt=h/m;
		 t=x;
		 for(k=0;k<m;k++)
		 {
			 gillf(t,y,(double *)d,n);
			 for(ii=0;ii<4;ii++)
			 {
				 for(j=0;j<n;j++)
				 {
					 d[j]=d[j]*hh;
					 for(j=0;j<n;j++)
					 {
						 r=a[ii]*(d[j]-b[ii]*g[j]);//是否有问题？
						 *(y+j)=*(y+j)+r;
						 s=g[j]+3.0*r;
						 g[j]=s-c[ii]*d[j];
					 }
					 t0=t+e[ii]*hh;
					 gillf(t0,y,(double *)d,n);
				 }
				 t=t+dt;
			 }
			 p=0.0;
			 for(j=0;j<n;j++)
			 {
				 qq=fabs(*(y+j)-v[j]);
				 if(qq>p)
					 p=qq;
			 }
			 hh=hh/2.0;
			 m=m+m;
		 }
		 for(j=0;j<n;j++)
			 q[j]=g[j];
		 return;
	}
}
*/

void equat_gill(double t,double h,double eps,double *y,double *q,int n)
{
	extern void gillf(double t,double *y,double *d,int);
	int i,j,k,m,ii;
	double x,p,hh,r,s,t0,dt,qq;
	double d[20],u[20],v[20],g[20];
	//求y[i]的一些参数
	double a[4]={0.5,0.29289321881,1.7071067812,0.166666667};
	double b[4]={2.0,1.0,1.0,2.0};//y[i]中q的系数
	double c[4];//q[i]中k的系数(1)
	double e[4]={0.5,0.5,1.0,1.0};//q[i]中k的系数(2)
	for(i=0;i<3;i++)
		c[i]=a[i];
	c[3]=0.5;
	x=t;
	p=1.0+eps;
	hh=h;
	m=1;
	for(j=0;j<n;j++)
		u[j]=*(y+j);//保存初始值
	while(p>=eps)
	{
		for(j=0;j<n;j++)
		{
			v[j]=*(y+j);
			*(y+j)=u[j];
			g[j]=q[j];
		}
		dt=h/m;
		t=x;
		for(k=0;k<m;k++)
		{
			gillf(t,y,(double *)d,n);
			for(ii=0;ii<4;ii++)
			{
				for(j=0;j<n;j++)
					d[j]=d[j]*hh;
				for(j=0;j<n;j++)
				{
					r=a[ii]*(d[j]-b[ii]*g[j]);
					*(y+j)=*(y+j)+r;
					s=g[j]+3.0*r;
					g[j]=s-c[ii]*d[j];
				}
				t0=t+e[ii]*hh;
				gillf(t0,y,(double *)d,n);
			}
			t=t+dt;
		}
		p=0.0;
		for(j=0;j<n;j++)
		{
			qq=fabs(*(y+j)-v[j]);
			if(qq>p)
				p=qq;
		}
		hh=hh/2;
		m=m+m;
	}
	for(j=0;j<n;j++)
		q[j]=g[j];
	return ;
}

int equate_qr(double *a,double *xr,double *xi,double eps,int n,int jt)
{
	int i,j;
	double q[100];//方程的最大阶数限制在10
    for(j=0;j<n;j++)
		q[j]=-*(a+n-j-1)/(*(a+n));
	for(j=n;j<n*n;j++)
		q[j]=0.0;
	for(i=0;i<n-1;i++)
		q[(i+1)*n+i]=1.0;
	i=matrix_hqr((double *)q,xr,xi,eps,n,jt);
	return i;
}