C++一元三次方程求解算法 - 小众知识

C++一元三次方程求解算法

2013年01月27日 14:18:05 苏内容
  标签: C++/算法
阅读:7328

    x = z - b/3a,
代入可得
    z^3 + pz + q = 0,
对上面化简的方程,有求根公式:

x = (-q / 2 + (q ^ 2 / 4 + p ^ 3 / 27) ^ (1 / 2)) ^ (1 / 3) + (-q / 2 - (q ^ 2 / 4 + p ^ 3 / 27) ^ (1 / 2)) ^ (1 / 3),

x = w (-q / 2 + (q ^ 2 / 4 + p ^ 3 / 27) ^ (1 / 2)) ^ (1 / 3) + w ^ 2 (-q / 2 - (q ^ 2 / 4 + p ^ 3 / 27) ^ (1 / 2)) ^ (1 / 3),

x = w ^ 2 (-q / 2 + (q ^ 2 / 4 + p ^ 3 / 27) ^ (1 / 2)) ^ (1 / 3) + w (-q / 2 - (q ^ 2 / 4 + p ^ 3 / 27) ^ (1 / 2)) ^ (1 / 3).

其中w是三次本原单位根
w = ( i√3 - 1 ) / 2.


4.用牛顿迭代法比较好(特别是解高阶方程的时候),比如用牛顿迭代法求方程 
  2x3-4x2+3x-6=0   在1.5附近根。 代码如下:  
  #include   <math.h>   
  main()   
  {   
    float   x,x0,f,f1;   /*   f指左边的函数,f1指f的一阶倒数*/     
      x=1.5;   
      do   
    {       x0=x;   
              f=((2*x0-4)*x0+3)*x0-6;   
              f1     =(6*x0-8)*x0+3;   
              x=x0-f/f1;   
      }   
      while   (fabs(x-x0)>=1e-5);   
    printf("方程的根=%6.2f\n",x);   
  }
用二分法求方程在(-10,10)之间的根2X(3)-4X(2)+3X-6=0
#include<stdio.h>   
#include<math.h>   
  main()   
  {   
  float   x0,x1,x2,fx0,fx1,fx2;   
  do   
  {   
    printf("Please   input   x1   x2   value:\n");   
    scanf("%f,%f",&x1,&x2);   
    fx1=x1*((2*x1-4)*x1+3)-6;   
    fx2=x2*((2*x2-4)*x2+3)-6;   
  }   
  while(fx1*fx2>0);   
  do   
  {   
    x0=(x1+x2)/2;   
    fx0=x0*((2*x0-4)*x0+3)-6;   
    if((fx0*fx1)<0)   
    {   
      x2=x0;   
      fx2=fx0;   
    }   
    else    
    {   
      x1=x0;   
      fx1=fx0;   
    }   
  }   
  while(fabs(fx0)>=1e-5);   
  printf("The   fachengdegen   is   =%6.2f\n",x0);


# include <iostream.h>
# include <stdio.h>
# include <math.h>

//////////////////函数--用卡丹公式解一元三次方程/////////////////
int fun(double a,double b,double c,double d,
   double &real_y1,double &real_y2,double &real_y3,
   double &imag_y1,double &imag_y2,double &imag_y3)
{
 double p,q,r,u,v,g,h,fai;
 if (a ==0)
 return -1;
 p=(3.0*a*c-b*b)/(3*a*a);
 q=(2.0*pow(b,3.0)-9*a*b*c+27.0*a*a*d)/(27.0*pow(a,3.0));
 r=b/(3.0*a);
 h=pow(q/2.0,2.0)+pow(p/3.0,3.0);
 g=sqrt(h);
 if(h>=0)
 {
  if(-q/2.0+g<0)
  u=-pow(fabs(-q/2.0+g),1.0/3.0);
  else
  u=pow((-q/2.0+g),1.0/3.0);
  if(-q/2.0-g<0)
  v=-pow(fabs(-q/2.0-g),1.0/3.0);
  else
  v=-pow((-q/2.0-g),1.0/3.0);
  if(h==0)
  {
   real_y1=u+v-r;            imag_y1=0;
   real_y2=-(u+v)/2-r;       imag_y2=0;
   real_y3=-(u+v)/2-r;       imag_y3=0;
  }
  else
  { 
   real_y1=u+v-r;       imag_y1=0;
   real_y2=-(u+v)/2;    imag_y2=sqrt(3.0)*(u-v)/2;
   real_y3=-(u+v)/2;    imag_y3=-sqrt(3.0)*(u-v)/2;
  }
 }
 else
 {
  fai=acos((-q/2)/(sqrt(pow(fabs(p),3)/27)));
  real_y1=2*sqrt(fabs(p)/3.0)*cos(fai/3.0)-r;
  real_y2=-2*sqrt(fabs(p)/3.0)*cos((fai+3.1415926)/3.0)-r;
  real_y3=-2*sqrt(fabs(p)/3.0)*cos((fai-3.1415926)/3.0)-r;
  imag_y1=0;   imag_y2=0;    imag_y3=0;
 }
 return 0;
}
//////////////////////////////////主函数////////////////////////////////
 void main()
 {
  double a,b,c,d;
  double real_x1,real_x2,real_x3;
  double imag_x1,imag_x2,imag_x3;

 /* cout<<"请输入方程的系数a,b,c,d:"<<"\n"<<endl;
  cout<<"系数a=";
  cin>>a;
  cout<<endl;
  cout<<"系数b=";
  cin>>b;
  cout<<endl;
  cout<<"系数c=";
  cin>>c;
  cout<<endl;
  cout<<"系数d=";
  cin>>d;
  cout<<endl;*/


  //cubic 1
a = 18209.7;
b = 607922;
c = -1.06899e+09;
d = 1.63462e+08;
  if (fun(a,b,c,d,real_x1,real_x2,real_x3,imag_x1,imag_x2,imag_x3)!=-1)
  {
  cout<<"  "<<"   "<<"根的实部"<<"      "<<"根的虚部"<<"\n"<<endl;
  printf("x1   %.5f      %.5f\n\n",real_x1,imag_x1);
  printf("x2   %.5f      %.5f\n\n",real_x2,imag_x2);
  printf("x3   %.5f      %.5f\n\n",real_x3,imag_x3);
  }
  else
  {
  cout << "not cubic equation" << endl;
  }


  //cubic 2
a = 22116.5;
b = 266288;
c = -5.32973e+08;
d = 1.76198e+08;
   if (fun(a,b,c,d,real_x1,real_x2,real_x3,imag_x1,imag_x2,imag_x3)!=-1)
  {
  cout<<"  "<<"   "<<"根的实部"<<"      "<<"根的虚部"<<"\n"<<endl;
  printf("x1   %.5f      %.5f\n\n",real_x1,imag_x1);
  printf("x2   %.5f      %.5f\n\n",real_x2,imag_x2);
  printf("x3   %.5f      %.5f\n\n",real_x3,imag_x3);
  }
  else
  {
  cout << "not cubic equation" << endl;
  }


  //cubic 3
a = 26968.2;
b = 315403;
c = -5.9352e+08;
d = 2.0085e+08;
   if (fun(a,b,c,d,real_x1,real_x2,real_x3,imag_x1,imag_x2,imag_x3)!=-1)
  {
  cout<<"  "<<"   "<<"根的实部"<<"      "<<"根的虚部"<<"\n"<<endl;
  printf("x1   %.5f      %.5f\n\n",real_x1,imag_x1);
  printf("x2   %.5f      %.5f\n\n",real_x2,imag_x2);
  printf("x3   %.5f      %.5f\n\n",real_x3,imag_x3);
  }
  else
  {
  cout << "not cubic equation" << endl;
  }



 }


3次方程总是有实根,所以可以先使用牛顿迭代法计算出任意一个实根
比如
ax^3+bx^2+cx+d=0
先计算出实根x1
然后通过多项式除法计算
(ax^3+bx^2+cx+d)/(x-x1)=ax^2+ex+f
然后对方程
ax^2+ex+f=0使用根的判别式判断是否有实根就可以了。

使用卡丹公式只有理论上面的意义。
扩展阅读