C++求解一元三次方程的实根

最新推荐文章于 2024-06-28 08:48:13 发布

残影丶

最新推荐文章于 2024-06-28 08:48:13 发布

阅读量1w

点赞数 11

分类专栏： C++ matlab 文章标签： C++求解一元三次方程盛金公式求解一元三次方程

本文链接：https://blog.csdn.net/liumangmao1314/article/details/88796429

版权

C++ 同时被 2 个专栏收录

6 篇文章 2 订阅

订阅专栏

matlab

6 篇文章 4 订阅

订阅专栏

引言

一元三次方程的解法有多种，百度百科也有介绍，卡尔丹公式和盛金公式。我这里主要是用著名的盛金公式解一元三次方程，不得不说盛金公式更为直观性，与我们解一元二次方程很类似。下面就拿来直接用了。

盛金公式

$ax^{3}+bx^{2}+cx+d=0(a,b,c,d\in \mathbb{R},a\neq 0)$

重根判别式： $\left\{\begin{matrix} A=b^{2}-3ac\\ B=bc-9ac\\ C=c^{2}-3bd\\ \end{matrix}\right.$ ，总判别式为： $\triangle =B^{2}-4AC$ 。还有四个求根公式详见给的链接，这里就不一一写出了。

代码实现
C++

从上可以看出，转成C++/matlab很容易实现，但是注意C++判断一个数（float）等于0不能用==0，而是用<1e-6，所有有时候会有一些很小的误差。下面我给出我的实现函数ShengJin，调用即可。注意这里我只取出了实根。

void ShengJin(double a,double b,double c,double d,vector<double> &X123)
{
	/************************************************************************/
	/* 盛金公式求解三次方程的解 
	   德尔塔f=B^2-4AC
           这里只要了实根，虚根需要自己再整理下拿出来
	*/
	/************************************************************************/
	double A=b*b-3*a*c;
	double B=b*c-9*a*d;
	double C=c*c-3*b*d;
	double f=B*B-4*A*C;
	double i_value;
	double Y1,Y2;
	if (fabs(A)<1e-6 && fabs(B)<1e-6)//公式1
	{
		X123.push_back(-b/(3*a));
		X123.push_back(-b/(3*a));
		X123.push_back(-b/(3*a));
	}
	else if (fabs(f)<1e-6)   //公式3
	{
		double K=B/A;
		X123.push_back(-b/a+K);
		X123.push_back(-K/2);
		X123.push_back(-K/2);
	}
	else if (f>1e-6)      //公式2
	{
		Y1=A*b+3*a*(-B+sqrt(f))/2;
		Y2=A*b+3*a*(-B-sqrt(f))/2;
		double Y1_value=(Y1/fabs(Y1))*pow((double)fabs(Y1),1.0/3);
		double Y2_value=(Y2/fabs(Y2))*pow((double)fabs(Y2),1.0/3);
		X123.push_back((-b-Y1_value-Y2_value)/(3*a));//虚根我不要
		//虚根还是看看吧，如果虚根的i小于0.1，则判定为方程的一根吧。。。
		i_value=sqrt(3.0)/2*(Y1_value-Y2_value)/(3*a);
		if (fabs(i_value)<1e-1)
		{
			X123.push_back((-b+0.5*(Y1_value+Y2_value))/(3*a));
		}
	}
	else if (f<-1e-6)   //公式4
	{
		double T=(2*A*b-3*a*B)/(2*A*sqrt(A));
		double S=acos(T);
		X123.push_back((-b-2*sqrt(A)*cos(S/3))/(3*a));
		X123.push_back((-b+sqrt(A)*(cos(S/3)+sqrt(3.0)*sin(S/3)))/(3*a));
		X123.push_back((-b+sqrt(A)*(cos(S/3)-sqrt(3.0)*sin(S/3)))/(3*a));
	}
}

这里参照博客matlab实现一元三次方程来验证代码。

matlab

matlab有自带的函数多项式求解函数：roots 求解一元三次方程。

x=roots([a, b, c, d]);
for k=1:length(x)
    xx(k) = isreal(x(k));%求实根
end
X_real=x(xx); %X_real都是实根

当然，matlab也有大佬实现了盛金公式求解一元三次方程，可以参考这篇博客，写的很详细：盛金公式matlab版。

总结

在编写求解代码时，也看过和实践过别人的求解代码，很多都是用基于一次二次导数的二分法来求解方程（我还有一份，需要的可以私信我），较为复杂且费时。因此参照盛金公式进行了编写C++代码，网上肯定也有同类代码，思想都一样的，很多时候我们为了求解一元三次方程并不是实验，而是为了达到我们的研究目的，比如我这只要实根。