c++:求解二元二次方程组（解析解）

这个问题从产生到解决，困扰了LZ快一周（当然因为中间也有其他事情要做），本来从LZ的角度是想找第三方库使用，不要说matlab就可以啊o(╥﹏╥)o，确实matlab一个solve函数就解决问题了，但是问题是不能用matlab，所以第一方案排除。

用第三方库？google了半天，发现一个gsl库确实可以求解二元二次方程组，但是也有问题，1）只能求解出一个解，而且无法求解复数解（源代码定义中没有复数定义，除非自己扒源程序，自己写）；2）使用的是优化算法，例如牛顿迭代等，所求出的是近似解，非解析解，也就是说求得的结果并不精确。而且初始值的设定会影响迭代结果，容易进入局部最优。3）复数解无法收敛。

PS：花钱的库还是算了，毕竟也不是特别难解决的问题。。。

只剩一个方案：自己动手，丰衣足食O(∩_∩)O哈哈~

列出基本的问题：
a1x2+b1xy+c1y2+d1x+e1y+f1=0a_1x^2+b_1xy+c_1y^2+d_1x+e1y+f_1=0a1x2+b1xy+c1y2+d1x+e1y+f1=0

a2x2+b2xy+c2y2+d2x+e2y+f2=0a_2x^2+b_2xy+c_2y^2+d_2x+e2y+f_2=0a2x2+b2xy+c2y2+d2x+e2y+f2=0

求上述方程的解（包含复数），也就是两个二次曲线的交点。

step 1: 利用两式相减，用x来表示y
得到类似如下的方程：

(Bx+D)y=−Ax2−Cx−E(Bx+D)y = -Ax^2-Cx-E(Bx+D)y=−Ax2−Cx−E

其中的A,B,C,D,E都是用a1−f1,a2−f2a1-f1 , a2-f2a1−f1,a2−f2 所表示的。

step 2:带入原方程得到关于xxx的一元四次方程

M4x4+M3x3+M2x2+M1x+M0=0M_4x^4 + M_3x^3+M_2x^2 + M_1x +M_0 = 0M4x4+M3x3+M2x2+M1x+M0=0

将四次项归一化

x4+bx3+cx2+dx+e=0x^4+bx^3+cx^2+d^x+e=0x4+bx3+cx2+dx+e=0

step 3: 利用费拉里法求一元四次方程

具体对费拉里法求解一元四次方程的推导：https://baike.baidu.com/item/%E4%B8%80%E5%85%83%E5%9B%9B%E6%AC%A1%E6%96%B9%E7%A8%8B

这里就不在解释了，因为公式很多，手打太浪费时间了。

讲了这么多，贴出源代码：（如果有小伙伴有更好的方法可以和LZ私聊哦）

#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <complex>
// #include <cmath>
#include <math.h>
#include <cfloat>void quadratic2quartic(Eigen::Matrix<double, 6, 1> & _par1, Eigen::Matrix<double, 6, 1> & _par2,Eigen::Matrix<double, 5, 1> & _y2x_par,Eigen::Matrix<std::complex<double>, 5, 1> &_quartic_par)
{double A, B, C, D, E;double a1, b1, c1, d1, e1, f1, a2, b2, c2, d2, e2, f2;a1 = _par1[0];b1 = _par1[1];c1 = _par1[2];d1 = _par1[3];e1 = _par1[4];f1 = _par1[5];//   std::cout  << a1 << "   " << b1 <<  "   " <<  c1  << "   " << d1 <<  "   " << e1 <<  "   " << f1 << std::endl;a2 = _par2[0];b2 = _par2[1];c2 = _par2[2];d2 = _par2[3];e2 = _par2[4];f2 = _par2[5];//   std::cout  << a2 << "   " << b2 <<  "   " <<  c2  << "   " << d2 <<  "   " << e2 <<  "   " << f2 << std::endl;A = a1*c2 - c1*a2;B = b1*c2 - c1*b2;C = d1*c2 - c1*d2;D = e1*c2 - c1*e2;E = f1*c2 - c1*f2;_y2x_par << A, B, C, D, E;
//   std::cout << "A-E: \n";
//   std::cout << A << "   " << B << "   " << C <<"   " << D << "   " << E << std::endl; double M0, M1, M2, M3, M4;M4 = a1*B*B - b1*A*B + c1*A*A;M3 = 2*a1*B*D - b1*A*D - b1*B*C + 2*c1*A*C + d1*B*B - e1*A*B;M2 = a1*D*D - b1*C*D - b1*B*E + c1*C*C + 2*c1*A*E + 2*d1*B*D - e1*A*D - e1*B*C + f1*B*B;M1 = -b1*D*E + 2*c1*C*E + d1*D*D - e1*C*D - e1*B*E + 2*f1*B*D;M0 = c1*E*E - e1*D*E + f1*D*D;
//   std::cout << "m4-m0: \n";
//   std::cout << M4 << "   " << M3 << "   " << M2 << "   " << M1 << "   " << M0 << std::endl;
//   double b, c, d, e;b = M3 / M4;c = M2 / M4;d = M1 / M4;e = M0 / M4;
//   std::cout  << "b, c, d, e: \n";
//   std::cout  << " " << b << "   " << c << "   " << d<< "   " << e << "   \n"; _quartic_par << 1.0, b, c, d, e;}//extract a root
void sqrtn(const std::complex<double> &_input, double n,std::complex<double> & _out)
{
//   double r = sqrt(_input.real()*_input.real() + _input.imag()*_input.imag());
//   std::cout << _input << std::endl;double r = hypot(_input.real(), _input.imag());//   std::cout << "the norm is " << r << std::endl;if(r > 0.0){double a = atan2(_input.imag(), _input.real());
//     std::cout << "the a is: " << a << std::endl;//     double a = arg(_input);n = 1 / n;r = pow(r, n);a *= n;_out.real() = r * cos(a) ;_out.imag() = r * sin(a);}else{_out.real() = 0.0;_out.imag() = 0.0;}
//   std::cout << "out :" << _out << std::endl;
}void Ferrari(Eigen::Matrix<std::complex<double>, 5, 1> & _quartic_par, Eigen::Matrix<std::complex<double>, 4, 1> & _x)
{std::complex<double> a = _quartic_par[0];std::complex<double> b = _quartic_par[1];std::complex<double> c = _quartic_par[2];std::complex<double> d = _quartic_par[3];std::complex<double> e = _quartic_par[4];std::complex<double> P = (c*c + 12.0*e - 3.0*b*d) / 9.0;std::complex<double> Q = (27.0*d*d + 2.0*c*c*c + 27.0*b*b*e - 72.0*c*e - 9.0*b*c*d) / 54.0;
//   std::cout << "P: " << P << "   \nQ: " << Q << std::endl;std::complex<double> D, u, v;
//   D = cabs(sqrt(Q*Q - P*P*P));sqrtn(Q*Q - P*P*P, 2.0, D);
//   std::cout << "D: " << D << std::endl;u = Q + D;v = Q - D;
//   std::cout << "u: " << u << "  \nv: " << v << std::endl;if(v.real()*v.real() + v.imag()*v.imag() > u.real()*u.real() + u.imag()*u.imag()){sqrtn(v, 3.0, u);}else{sqrtn(u, 3.0, u);}
//   std::cout <<"u: " << u << std::endl;std::complex<double> y;if(u.real()*u.real() + u.imag()*u.imag() > 0.0){v = P / u;std::complex<double> o1(-0.5,+0.86602540378443864676372317075294);std::complex<double> o2(-0.5,-0.86602540378443864676372317075294);std::complex<double> &yMax = _x[0];double m2 = 0.0;double m2Max = 0.0;int iMax = -1;for(int i = 0; i < 3; ++i){y = u + v + c / 3.0;u *= o1;v *= o2;a = b*b + 4.0*(y-c);m2 = a.real()*a.real() + a.imag()*a.imag();if(0==i || m2Max < m2){m2Max = m2;yMax = y;iMax = i;}}y = yMax;}else{//cubic equationy = c / 3.0;  }std::complex<double> m;sqrtn(b*b + 4.0*(y-c), 2.0, m);if(m.real()*m.real() + m.imag()*m.imag() >= DBL_MIN){std::complex<double> n = (b*y - 2.0*d) / m;sqrtn((b+m)*(b+m) - 8.0*(y+n), 2.0, a);_x[0] = (-(b+m) + a) / 4.0;
//     std::cout << "_x[0]: " << _x[0] << std::endl;_x[2] = (-(b+m) - a) / 4.0;
//     std::cout << "_x[1]: " << _x[1] << std::endl;sqrtn((b-m)*(b-m) - 8.0*(y-n), 2.0, a);_x[1] = (-(b-m) + a) / 4.0;_x[3] = (-(b-m) -a) / 4.0;
//     std::cout << "_x[2]: " << _x[2] << std::endl;
//     std::cout << "_x[3]: " << _x[3] << std::endl;}else{sqrtn(b*b - 8.0*y, 2.0, a);_x[0] = _x[1] = (-b + a) / 4.0;_x[2] = _x[3] = (-b - a) / 4.0;}//   for(int i = 0; i <4; i++)
//   {
//     std::cout << _x[i] << std::endl;
//   }
}void compute_y(Eigen::Matrix<std::complex<double>, 4, 1> &_x, Eigen::Matrix<double, 5, 1> &_y2x_par, Eigen::Matrix<std::complex<double>, 4, 1> &_y)
{for(int i = 0; i < 4; i++){_y[i] = (-_y2x_par[0]*_x[i]*_x[i] - _y2x_par[2]*_x[i] - _y2x_par[4]) / (_y2x_par[1]*_x[i] + _y2x_par[3]);
//     std::cout << "_y[" << i << "]" << _y[i] << std::endl;}
}int main(int argc, char **argv) {Eigen::Matrix<double, 6, 1> ic1, ic2;ic1 <<414181, -12635.8, 398758, -9.53484e+08, -4.16931e+08, 5.00382e+11;ic2 << 923938, -41764.7, 932015, -2.10858e+09, -9.79279e+08, 6.2537e+11;Eigen::Matrix<std::complex< double >, 5, 1> quartic_par;Eigen::Matrix<double, 5, 1> y2x_par;
//   quadratic2quartic(ic1, ic2, quartic_par);quadratic2quartic(ic1, ic2, y2x_par, quartic_par);Eigen::Matrix<std::complex< double >, 4, 1> x, y;Ferrari(quartic_par, x);compute_y(x, y2x_par, y);std::cout <<"the result :" << std::endl;for(int i = 0; i < 4; i++){std::cout << "x: " << x[i] <<  "\ty: " << y[i] << std::endl;}//     std::cout << "Hello, world!" << std::endl;return 0;
}

使用了Eigen处理矩阵会比较方便。

贴出最后运行结果：

使用matlab对上述算法进行验证：

最后运行结果如下所示：

好，最后检验完毕。最后感慨一下还是matlab功能强大，如果要使用c++还是要对具体算法原理了解，才能编写代码，对编程能力、算法理解要求更高，所以小伙伴们自行决定需要什么工具哦O(∩_∩)O哈哈~解决一个问题好开心(｡◕ˇ∀ˇ◕)

c++:求解二元二次方程组（解析解）相关推荐

c++编程求解二元二次方程组_c++ 牛顿迭代法求解多元多次方程组
//经典牛顿迭代法C++实现 #include #include #define N 3 // 非线性方程组中方程个数.未知量个数 #define Epsilon 0.000000001 // 差向量 ...
c++编程求解二元二次方程组_一道俄罗斯高难度解方程组题，错误率达99%+，中国学霸：确实很难...
同学们好,今天老师为大家分享一道俄罗斯奥林匹克竞赛题之高难度解方程组试题.提到解方程组,同学们第一反应就是想到之前我们学过的一些解方程组的方法与技巧,例如:加减消元法.代入消元法.换元法等.也就是说, ...
c++编程求解二元二次方程组_C++编程风格约定
序 C++用法很多,包容性也比较强.一个C++的工程可能包含了各种各样没见过的用法.本篇内容主要是参照谷歌C++标准规范,结合自身实际工作及经验,整理一份适合平时C++开发的规则,规范自身C++编程 ...
python计算二元二次方程组_2020高中数学初高中衔接读本专题4.1简单的二次方程组的解法精讲深剖学案202020211109...
在初中我们已经学习了一元一次方程.一元二次方程及二元一次方程组的解法,掌握了用消元法解二元一次方程组．高中学习圆锥曲线时,需要用到二元二次方程组的解法．因此,本讲讲介绍简单的二元二次方程组的解法． ...
【原创】《矩阵的史诗级玩法》连载三十二：用矩阵法解二元二次方程组的一般式
现在我们给出一个方程组,然后尝试用矩阵来求解. 在连载十六中,我们给出了曲线类型的判断法则: Δ<0时,方程为椭圆(包括正圆) Δ>0时,方程为双曲线 Δ=0时,方程为抛物线其中Δ=B^ ...
python求解三元一次方程_北师大版八上数学5.2 求解二元一次方程组知识点微课精讲...
知识点总结代入消元法代入消元法的实质是将二元一次方程组中的某一个方程进行未知数的分离,即将该方程进行变换,完整分离出一个独立的未知数,而这个未知数将用含有另一个未知数的式子来表示.设某二元一次方程 ...
二元二次方程例题_二元二次方程组例题_相关文章专题_写写帮文库
时间:2019-05-15 01:34:12 作者:admin 初三代数教案第十二章:一元二次方程第20课时:由一个二元一次方程和一个二元二次方程组成的方程组(一) 教学目标: 1.使学生了解二 ...
C语言求解一元二次方程组的代码
#include <stdio.h> #include <math.h>int main() {double a, b, c, deta, x1, x2, p, q;scanf ...
matlab使用solve求解二元二次方程组
网上有些代码存在问题,这里只做了原理上的更正.也希望能各位码农能完全.准确地提供代码. 求解二元二次方程组代码: clc clear close all syms x y; f_1 = sym(x^2 ...

c++:求解二元二次方程组（解析解）

c++:求解二元二次方程组（解析解）相关推荐

最新文章

热门文章