计算多项式非线性方程f(x)=0的求根问题 .cpp

来自「计算多项式非线性方程f(x)=0的求根问题,本人数值分析课程设计做的一个程序,供」· C++ 代码 · 共 489 行

/***************************实验题目:计算非线性方程f(x)=0的求根问题***********************


一:设计内容和期望
     本程序把函数f(x)在X(k) 完成二次泰勒展开,并得到一个新的迭代公式计算X(k+1),以此来计算非线性方程
     f(x)=0的根,从而改进牛顿法的计算公式,使计算精度从平方收敛达到3阶收敛(估计,尚待证明)


二:遇到的问题
    1:对于任意f(x)的输入;
	2:f(x)导数的实现;
	3:误差精度的控制;


三:根据以上问题,提出下列解决方案
    1:采用数据结构的相关知识解决输入问题,(由于受到算法限制先完成对f(x)为多项式函数的形式,其他
	    例如包含sin,cos,exp的形式有待解决);
    2:编写通用多项式求导和求值问题,以便完成新的迭代公式;
	3:控制x(k+1)与x(k)的误差(精度可手动输入控制),保证迭代数列收敛,以便达到精度为3阶收敛;
	

四:实验效果预期
    1:可手动输入多项式形式的非线性方程;
	2:可手动输入精度,迭代初值;
	3:通过与牛顿法进行对比分析,从而从实验中证明精度有所提高;
	  (例如同样的多项式非线性方程,在同样的精度下,新方法比牛顿法迭代次数少即可证明精度提高)*/






#include "stdio.h"
#include "iostream.h"
#include "stdlib.h"
#include "math.h"
#include "iomanip.h"
/**********************************
*  	常量定义
*********************************/


#define TRUE 1
#define FALSE 0
#define OK 1
#define ERROR 0


/**********************************
*  	类型定义
*********************************/

typedef int Status;  /* 函数返回类型   */

/*  链表结点的数据元素的数据类型定义 */

typedef struct{
  float coef;
  int expn;
} Term,LkLtElemtype;

/*  链表结点的数据类型定义 */

typedef struct LNode{
  LkLtElemtype data;         /* 说明结点的数据元素类型为上面定义的结构体类型 */
  struct LNode *next;
} LNode,*Link, *Position;


/* 链表的数据类型定义 */

typedef struct{
  Link head,tail;   /*  head指向链表的头结点，tail指向最后一个结点 */
  int len;         /* 链表的长度*/
}LinkList;

/* 定义一个多项式类型为链表类型  */

typedef LinkList Polynomial;


/********************************** 
*   对链表的操作实现
*********************************/


/* 分配由*p指向的值为e的结点，并返回OK, 若分配失败返回ERROR  */

Status MakeNode(Link *p,LkLtElemtype e){

  *p=(Link) malloc(sizeof(LNode));
  if(!(*p)) return ERROR;
  (*p)->data=e;
  (*p)->next=NULL;
  return OK;
}


/* 释放*p所指向的结点  */

void FreeNode(Link *p){
  if(!(*p)){
    free(*p);
    *p=NULL;
  }
}


/*  构造一个空的线性链表,成功返回OK,失败返回ERROR  */

Status InitList(LinkList * l){
  if(!l) return ERROR;
  l->head=(Link)malloc(sizeof(LNode));
  if(!(l->head)) return ERROR;
  l->head->next=NULL;
  l->tail=NULL;
  l->len=0;
  return OK;
}


/* 销毁线性表链表*/

Status DestroyList(LinkList * l){
  Link q=NULL;
  Link p=NULL;
  if(!l) return ERROR;
  p=l->head->next;
  while(!p){
    q=p->next;
    free(p);
    p=q;
  }
  free(l->head);
  return OK;
}

/* 已知h指向线性链表的头结点，将s所指向结点插入在第一个结点之前 */

Status InsFirst(Link head,Link s){
  if ((!head ) || (!s)) return ERROR;
  s->next=head->next;
  head->next=s;
  return OK;
}

/*  已知*p指向线性链表中的一个结点，用e更新*p所指结点中数据元素的值  */

Status SetCurElem(Link *p,LkLtElemtype e){
  if (!(*p)) return ERROR;
  (*p)->data=e;
  return OK;
}


/*  已知p指向线性链表中的一个结点，返回p所指结点中数据元素的值  */

LkLtElemtype GetCurElem(Link p){
  return p->data;
}


/* 返回线性链表中头结点的位置 */

Position GetHead(LinkList l){
   return l.head;
}

/*************************************************************
*   对多项式链表所特有的操作实现(见教材P42-43)
*************************************************************/

/*  依a的指数值小于、等于、大于b的指数值，分别返回-1,0,1 */


int cmp(Term a,Term b){
  if(a.expn<b.expn) return -1;
  if(a.expn==b.expn) return 0;
  if(a.expn>b.expn) return 1;
  return OK;
}


/* 若有序链表中存在与e满足判定函数compare()取值为0的元素，则*q批示第一个值为e的结点的位置，并返回TRUE；
*  否则*q批示第一个与e满足判定函数compare()取值>0的元素的前驱的位置，并返回FALSE
*/
Status LocateElem2(LinkList l,LkLtElemtype e,Position *q,int (*compare)(LkLtElemtype,LkLtElemtype)){
  Link p;
  Link pp;
  p=l.head;
  do{pp=p,p=p->next;}
   while(p&&(*compare)(p->data,e)<0);

  if  (!p||(*compare)(p->data,e)>0)
 {*q=pp;return FALSE;}
else  {*q=p;return TRUE;}
}

/* 输入m项的系数和指数，建立表示一元多项式的按指数有序的链表 */

void CreatePolyn(Polynomial *p,int m)
{
  Term *t;
  int i;
  Link q;
  Link s;
  Link h;
  float f;
  InitList(p);
 /* printf("Ployn length:%d\n",(*p).len);*/
  h=GetHead(*p);
 /* printf("head:%o",h);*/
  t=(Term *)malloc(sizeof(Term));
  t->coef=0.0f;
  t->expn=-1;
  SetCurElem(&h,*t);
/* printf("head:expn:%d\n",(h->data).expn);*/
  printf("请输入一元多项式f(x):(格式:系数 指数,顺序可任意)\n");
  for(i=1;i<=m;++i)
  {
    //printf("please input the %d coef:",i);
    scanf("%f",&f);
    t->coef=f;
    //printf("please input the %d expn:",i);
    scanf("%d",&(t->expn));
    if(!LocateElem2(*p,*t,&q,(*cmp))){
      if(MakeNode(&s,*t)) InsFirst(q,s);
    }
  }
}



/* 显示一个用链表表示的多项式 */


void PrintPolyn(Polynomial p){
  Link q;
  int i=1;
  float coef;
  int expn;
  q=p.head->next;
  while(q)
  {
      coef=(q->data).coef;
      expn=(q->data).expn;//cout<<coef<<"x^"<<expn<<endl;
      if(i%10==0) printf("\n");
      else
	  {
		  if(expn!=0)
		{
          if(coef>0)
		  {
            if(i==1)  printf("%3.2fx^%d",coef,expn);
            else      printf("+%3.2fx^%d",coef,expn);
          }
          else printf("%3.2fx^%d",coef,expn);
       }
       else
	   {
         if(coef>0)
		 {
            if(i==1)  printf("%3.2f",coef);
            else      printf("+%3.2f",coef);
          }
          else printf("%3.2f",coef);
		 }
	}
      
      q=q->next;
      i++;
  }
  printf("\n");
 
}


/*多项式求导函数*/

void qiudao(Polynomial *p,Polynomial pa,int m)
{
    Link q1;
	q1=pa.head->next;
	Term *t;
  //int i;
  Link q;
  Link s;
  Link h;
  
  InitList(p);
   h=GetHead(*p);
 /* printf("head:%o",h);*/
  t=(Term *)malloc(sizeof(Term));
  t->coef=0.0f;
  t->expn=-1;
  SetCurElem(&h,*t);
  while (q1)
  {
    t->coef=q1->data.coef*q1->data.expn;
    t->expn=q1->data.expn-1;
	if(t->coef!=0)
	{
       if(!LocateElem2(*p,*t,&q,(*cmp)))
	   {
         if(MakeNode(&s,*t)) InsFirst(q,s);
	   }
	}
	
	
	q1=q1->next;
  }
}

/*多项式求值*/

double qiuzhi(Polynomial p,double x)
{
	Link q;double sum=0.0;
    q=p.head->next;
	while(q)
	{
		sum=sum+q->data.coef*pow(x,q->data.expn);
		q=q->next;
	}
	return sum;
}

/*初值条件的判断*/

int czpanduan(double *x,Polynomial pa,Polynomial pb,Polynomial pc)
{
	double f,f1,f2,sq;
        f=qiuzhi(pa,*x);f1=qiuzhi(pb,*x);f2=qiuzhi(pc,*x);
        sq=f1*f1-2*f2*f;	//cout<<"sq="<<sq<<"  "<<"f1="<<f1<<"  "<<"f2="<<f2<<endl;
	if(sq>=0&&f2!=0&&f1!=0) return OK;
	else return ERROR;
	
}

/*利用新的迭代公式计算收敛数列x(k+1)*/

double xindegongshi(double x,Polynomial pa,Polynomial pb,Polynomial pc)
{
     double f,f1,f2;
	 double t,t1;    //cout<<"x="<<x;
	f=qiuzhi(pa,x);
    //printf("f(%2f)=%f\n",x,f);
      f1=qiuzhi(pb,x);
    //printf("f`(%2f)=%f\n",x,f1);
     f2=qiuzhi(pc,x);
    //printf("f``(%2f)=%f\n",x,f2);
    t=x-(f1+sqrt(f1*f1-2*f2*f))/f2;
	t1=x-(f1-sqrt(f1*f1-2*f2*f))/f2;//新的迭代公式

	//cout<<"t="<<t;cout<<"t="<<t1;
    
	if(fabs(t-x)>fabs(t1-x)) return t1;
	else return t;//通过与上一次的值比较,选择绝对值差小的作为下一次的迭代值,从而保证数列收敛
}




int xinfangfa(double x,double exp,Polynomial pa,Polynomial pb,Polynomial pc)
{ 
	
	double x1;int i=0;
	x1=xindegongshi(x,pa,pb,pc);
     if(!czpanduan(&x1,pa,pb,pc)) return ERROR;
	        cout<<"x["<<i++<<"]="<<setprecision(8)<<x<<endl;
			cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
     while (fabs(x1-x)>=exp)
	 {
		 x=x1;x1=xindegongshi(x,pa,pb,pc);
		 cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
	 }
	 return OK;
}

/*newton迭代公式*/

double gongshi(double x,Polynomial pa,Polynomial pb,Polynomial pc)
{
     double f,f1;
	 double t;   
	f=qiuzhi(pa,x);
    f1=qiuzhi(pb,x);
    t=x-f/f1;
    return t;
}


void newton(double x,double exp,Polynomial pa,Polynomial pb,Polynomial pc)
{ 
	
	double x1;int i=0;
	x1=gongshi(x,pa,pb,pc);
	        cout<<"x["<<i++<<"]="<<setprecision(8)<<x<<endl;
			cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
     while (fabs(x1-x)>=exp)
	 {
		 x=x1;x1=gongshi(x,pa,pb,pc);
		 cout<<"x["<<i++<<"]="<<setprecision(8)<<x1<<endl;
	 }
}





/*************************************************************
*               主函数
*************************************************************/

void main()
{

  Polynomial pa,pb,pc;/* 定义多项式链表pa,pa,pc分别保留原函数,一阶导函数,二阶导函数 */

  int ma; /*pa的项数*/

  double x;/*x作为初值*/
  
  double exp;/*定义精度*/
  
  /*输入多项式函数项数*/

  printf("请输入多项式函数项数:");
  scanf("%d",&ma);

  /*按指数有序构造多项式链表pa,pb,pc*/

  CreatePolyn(&pa,ma);qiudao(&pb,pa,ma);qiudao(&pc,pb,ma);

  /*显示多项式pa*/

  printf("你输入的多项式函数f(x)=");
  PrintPolyn(pa);
  
  /*显示多项式一阶和二阶导函数*/

  printf("f`(x)=");
  PrintPolyn(pb);
  printf("f``(x)=");PrintPolyn(pc);cout<<endl;

  /*输入迭代初值*/
 
  printf("现在计算f(x)=0的根,请输入一个迭代初值:");
  cin>>x;

  /*判断迭代初值是否满足要求*/
  while(1)
  {
   if(czpanduan(&x,pa,pb,pc)) break;
   else 
   {
	  printf("初值不和要求,请重新输入一个迭代初值:");
      cin>>x;
   }
  }
			  
  /*输入精度*/

  printf("请输入精度(例如1e-3,表示0.001)");
  cin>>exp;cout<<endl;

  /*利用新的迭代公式计算*/

  cout<<"利用新的迭代公式计算的收敛数列如下:"<<endl;
  while(1)
  {
	  if(xinfangfa(x,exp,pa,pb,pc)) break;
      printf("中间变量不和要求,请重新输入一个迭代初值,重新进行计算:");
      cin>>x;
  }
  cout<<endl;
     
 /*利用牛顿法的迭代公式计算*/

  cout<<"利用牛顿法的迭代公式计算的收敛数列如下:"<<endl;
  newton(x,exp,pa,pb,pc);cout<<endl;

  /*计算完毕后销毁多项式函数，释放空间*/
  DestroyList(&pa);DestroyList(&pb);DestroyList(&pc);
  
  
  cout<<"通过数据请分析该算法的优缺点,以便加以改进,谢谢!!"<<endl;
}