C语言--三次方程数值求解

三次方程计算器

一元三次方程的解析解法早在500年前就已经被意大利数学家攻克了：

我们可以看到，相比于二次方程，三次方程求根公式表达式要复杂的多，包含了一个二重根号，计算量庞大，而且大多数情况下会出现一对复数解，C语言不太好操作；在判别式小于0时，即使三个解都是实数，我们仍然需要对一个复数开立方，这样的一个公式显然是不适合人工计算和计算器编程的。

由于三次方程至少有一个实根，因此我们可以通过常规的数值方法找出这个实数解。本程序先用二分法估算方程零点的大致位置，再通过牛顿迭代法获得最大实数根的数值解，最后通过根与系数的关系：

将方程降为二次方程，通过二次方程的求根公式，获得该三次方程在复数范围内的全部解，这就避免了对复数开方，计算量也小得多。

本程序在理论上可以求解最大实根在[-64,+64]范围内的所有类型的三次方程。

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

double funeva(double a,double b,double c, double d, double x)
{
	double y=x*(x*(a*x+b)+c)+d;
	return y;
}
double diffeva(double a,double b,double c,double x)
{
	double y0=x*(3*a*x+2*b)+c;
	return y0;
}

double eqsolve(double a,double b,double c,double d)
{
   double e=-64,f=64;
   double x,m;
   int i=0; 
   while((f-e)>=0.25)
   {
   	 if(int(funeva(a,b,c,d,e))^int(funeva(a,b,c,d,(e+f)/2)<=0))
   	 {
   	 	f=(e+f)/2;
   	 	i++;
	 }
	 else
	 {
	 	e=(e+f)/2;
	 	i++;
	 }
   }
    x = (e+f)/2;
   do
   {
   	m=funeva(a,b,c,d,x)/diffeva(a,b,c,x);
   	x-= m;
   	i++;
   }
   while(fabs(m)>=10e-10);
   printf("i=%d\n",i);
   return x;
}

int main()
{
	double a,b,c,d,x1,x2,x3;
	double B,C,delta;
	double y;
	int i=0;
	while(i!=1)
	{
		printf("Please input the coefficient of cube equation:\n");
    	scanf("%lf,%lf,%lf,%lf",&a,&b,&c,&d);
	    x1=eqsolve(a,b,c,d);
    	printf("x1= %.8f\n",x1);
    	if(a!=0)
        {
         	B = b/a+x1;
         	C = -d/(a*x1);
    	    delta = B*B-4*C;
	        if(delta>=0)
            {
		         x2 = (-B+ sqrt(delta))/2.0;
		         x3 = (-B- sqrt(delta))/2.0;
		         printf("x2= %.8f\nx3= %.8f\n",x2,x3);
	        }
	        else
	       {
		         x2 = x3 = -B/2.0;
		    y = sqrt(-delta)/2.0;
	     	printf("x2= %.8f + %.8fi\n",x2,y);
		    printf("x3= %.8f - %.8fi\n",x2,y);
	        }
    	}
    	if(a==0)
    	{
    		x2 = -(c/b+x1);
    		printf("x2=%.8f\n",x2);
		}
	    printf("\n\nContinue 0(Yes)/1(No)?\n");
	    scanf("%d",&i);
	    printf("\n");
	    
	}
	system("pause");
	return 0;
}

运行结果：