一元四次方程求解-【附MATLAB代码】

目录

前言

求解方法

MATLAB验证

---

前言

        最近在研究机器人的干涉（碰撞）检测，遇到了一个问题，就是在求椭圆到原点的最短距离时，构建的方程是一个一元四次方程。无论是高中的初等数学，大学的高等数学，还是研究生的高等代数，都没有关于一元四次方程的求解方法，大多都是一元二次方程的求解。仔细一研究才知道为什么很少提及一元四次方程。具体解法如下：

求解方法

MATLAB验证

% output
% root：为方程的解（根），i 为解的个数
% input
% parameter 为方程的5个系数

function [root,i] = MYSolve4OrderEquaton(parameter)
a=parameter(2)/parameter(1);
b=parameter(3)/parameter(1);
c=parameter(4)/parameter(1);
d=parameter(5)/parameter(1);

a3=1;
b3=-b;
c3=(a*c-4*d);
d3=-(a^2*d-4*b*d+c^2);
parameter3=[a3,b3,c3,d3];
[root3,y3,i3] = Solve3OrderEquaton(parameter3);
i=0;
root=[];
for j=1:1
    if(a^2/4-b+root3(j)<0||root3(j)^2/4-d<0)
        break;
    end
    alpha=sqrt(a^2/4-b+root3(j));
    beta=sqrt(root3(j)^2/4-d);
    if(a*root3(j)/2-c>0)
        a21=1;
        b21=a/2-alpha;
        c21=root3(j)/2-beta;
        parameter21=[a21,b21,c21];
        [root21,y21,i21] = Solve2OrderEquaton(parameter21);
        a22=1;
        b22=a/2+alpha;
        c22=root3(j)/2+beta;
        parameter22=[a22,b22,c22];
        [root22,y22,i22] = Solve2OrderEquaton(parameter22);
    else
        a21=1;
        b21=a/2-alpha;
        c21=root3(j)/2+beta;
        parameter21=[a21,b21,c21];
        [root21,y21,i21] = Solve2OrderEquaton(parameter21);
        a22=1;
        b22=a/2+alpha;
        c22=root3(j)/2-beta;
        parameter22=[a22,b22,c22];
        [root22,y22,i22] = Solve2OrderEquaton(parameter22);
    end
    root4{j}=[root21,root22];
    i4{j}=[i21,i22];
    root=[root,root4{j}];
    i=i+i21+i22;
end
end
function [root,y,i] = Solve3OrderEquaton(parameter)
a=parameter(1);
b=parameter(2);
c=parameter(3);
d=parameter(4);
a_2=a*a;
a_3=a_2*a;
b_2=b*b;
b_3=b_2*b;
p=c/3/a-b_2/9/a_2;
q=d/2/a+b_3/27/a_3-b*c/6/a_2;
delta=q*q+p^3;
if(delta>0)
    i=1;
    root=nthroot(-q+sqrt(delta),3)+nthroot(-q-sqrt(delta),3)-b/3/a;
elseif(delta==0)
    i=2;
    root(1)=-2*nthroot(q,3)-b/3/a;
    root(2)=nthroot(q,3)-b/3/a;
else
    i=3;
    alpha=1/3*acos(-q*sqrt(-p)/p^2);
    root(1)=2*sqrt(-p)*cos(alpha)-b/3/a;
    root(2)=2*sqrt(-p)*cos(alpha+2/3*pi)-b/3/a;
    root(3)=2*sqrt(-p)*cos(alpha+4/3*pi)-b/3/a;
end
y=a*root.^3+b*root.^2+c*root+d;
end
function [root,y,i] = Solve2OrderEquaton(parameter)
a=parameter(1);
b=parameter(2);
c=parameter(3);
delta=b^2-4*a*c;
if(delta>0)
    i=2;
    root(1)=(-b+sqrt(delta))/2/a;
    root(2)=(-b-sqrt(delta))/2/a;
elseif(delta==0)
    i=1;
    root=-b/2/a;
else
    i=0;
    root=[];
end
y=a*root.^2+b*root+c;
end

测试代码

clc;
clear;
parameter = [1,1,1,1,0];
[root,i]=MYSolve4OrderEquaton(parameter);
[root,i]

结果验证