科学网 › 标签 › 一元四次方程

标签: 一元四次方程

相关帖子	版块	作者	回复/查看	最后发表

没有相关内容

相关日志

0058: 一元四次方程解法原理: cwhe10 2018-1-23 23:11; 一元四次方程解法原理作者：何成文 2018-1-23 1. 源代码（接口函数）： %---------------------------------------------------------------------- % 功能：一元四次方程解法 % 作者： Alex 2018-1-22 % 参考资料：百度百科：一元四次方程求根公式 % p= ; % 一元四次方程的解！ % xx=roots(p); %---------------------------------------------------------------------- function =Solve_Quartic_equation_Alex(a,b,c,d,e) deta1=c^2-3*b*d+12*a*e;deta2=2*c^3-9*b*c*d+27*a*d*d+27*b*b*e-72*a*c*e; wuyon=sqrt(-4*deta1^3+deta2^2)+deta2; deta=(deta1*2^(1/3))/(3*a*wuyon^(1/3))+ wuyon^(1/3)/(3*a*2^(1/3)); wuyon1=sqrt(b*b/(4*a*a) - 2*c/(3*a) + deta ); % 构造中间变量！ wuyon2=b*b/(2*a*a)- 4*c/(3*a) - deta ; wuyon3=(-b^3/a^3 + 4*b*c/a^2 -8*d/a )/(4*wuyon1); x1=-b/(4*a)-wuyon1/2-sqrt(wuyon2-wuyon3)/2; % 一元四次方程的四个根！ x2=-b/(4*a)-wuyon1/2+sqrt(wuyon2-wuyon3)/2; x3=-b/(4*a)+wuyon1/2-sqrt(wuyon2-wuyon3)/2; x4=-b/(4*a)+wuyon1/2+sqrt(wuyon2-wuyon3)/2; x= ; %---------------------------------------------------------------------- 2. 原理图片：; 个人分类: 科学研究|1 次阅读|0 个评论

又发现“费拉里”法使用的一个错误: 热度 1 Babituo 2012-6-3 16:08; 网上查到所谓的通用公式解法有“费拉里”法。基本想法是将一元四次方程的表达式配成两个完全平方表达式的等式，为了实现这个目的，需要在此之前解一个一元三次方程，得到配方的系数。解法如下：一般一元四次方程为： x^4+c1x^3+c2x^2+c3x+c4=0 移项后： x^4+c1x^3=-c2x^2-c3x-c4 两边同时加上(c1x/2)^2: x^4+c1x^3+(c1x/2)^2 = -c2x^2-c3x-c4+(c1x/2)^2 整理为： (x^2 + c1x/2)^2 = (c1^2/4-c2)x^2-c3x-c4 两边再同时加上 (x^2+c1x/2)y + (y/2)^2 (x^2 + c1x/2)^2 + (x^2+c1/2)y + (y/2)^2 = (c1^2/4-c2)x^2-c3x-c4 + (x^2+c1x/2)y + (y/2)^2 再整理： ((x^2+c1x/2）+ y/2)^2 = (c1^2/4-c2+y)x^2 + (c1y/2-c3)x +(y/2)^2-c4.....（1）左边为一个完全平方式，右边为一个x的二次多项式(系数中含有y), 为了使右边也成为一个完全平方式，必须使其判别式为零，即： DAT = ^2 - 4 * * = 0 于是，得到一个关于y的一元三次方程。 -y^3 + c2y^2 + (4c4-c1c3)y+c3^2+c4c1^2-4c4c2=0 只要解这个一元三次方程得到一个解y0，代入（1）右边，并令右边=0解出配方常量x0; 便可得到两边都为完全平方的方程，原一元四次方程被降为一元二次方程。可轻松求解。 ((x^2+c1x/2）+ y0/2)^2 =(x-x0)^2 即： (x^2+c1x/2）+ y0/2 = x-x0 (x^2+c1x/2）+ y0/2 = -x+x0 整理为： x^2+（c1/2-1）x + y0/2 + x0 = 0 x^2 +(c1/2+1）x + y0/2 - x0 = 0 解这两个一元二次方程的解，即为一元四次方程的解。以上是我使用费拉里法的过程，其中红色一行隐藏一个较难发现的错误。右边原来是： (c1^2/4-c2+y)x^2 + (c1y/2-c3)x +(y/2)^2-c4 当解出y=y0后代入，变为： (c1^2/4-c2+y0)x^2 + (c1y0/2-c3)x +(y0/2)^2-c4 令其为零，解出x=x0后，完全平方的表达式，不能约简平方项的系数，变为（x-x0）^2，应保留原有系数。令k=sqrt(c1^2/4-c2+y0) 完全平方式应写作：(kx-kx0)^2 这样，才能维持原1式表达式的左右相等。该行及之下应更正为：令k=sqrt(c1^2/4-c2+y0) ((x^2+c1x/2）+ y0/2)^2 =(kx-kx0)^2 即： (x^2+c1x/2）+ y0/2 = kx-kx0 (x^2+c1x/2）+ y0/2 = -kx+kx0 整理为： x^2+（c1/2-k）x + y0/2 + kx0 = 0 x^2 +(c1/2+k）x + y0/2 - kx0 = 0 解这两个一元二次方程的解，即为一元四次方程的解。用修改后的程序进行测试解方程：（x-1）*（x-2）*（x-3）*（x-4）= 0 即： x^4-10x^3+35x^2-50x+24 =0 得到： 4+0i 3+0i 2+0i 1+0i 测试通过！终于解决了一元四次方程的求解程序的问题，爽！; 个人分类: 信息探索|6498 次阅读|4 个评论

一元四次方程通用公式解法问题出在哪？: 热度 2 Babituo 2012-6-2 17:03; 在编写张学文老师的“时间胎”演示程序中，卡壳了。卡在求解根据屏幕坐标反投影求模型坐标时，要解一个一元四次方程。看起来并不复杂的一元四次方程，用迭代法很容易求解。但为了提高效率，想找个公式解法。网上查到所谓的通用公式解法有“费拉里”法。基本想法是将一元四次方程的表达式配成两个完全平方表达式的等式，为了实现这个目的，需要在此之前解一个一元三次方程，得到配方的系数。方法是不错，程序也编出来了，就是实际使用中怎么也得不到正确的解。解法如下：一般一元四次方程为： x^4+c1x^3+c2x^2+c3x+c4=0 移项后： x^4+c1x^3=-c2x^2-c3x-c4 两边同时加上(c1x/2)^2: x^4+c1x^3+(c1x/2)^2 = -c2x^2-c3x-c4+(c1x/2)^2 整理为： (x^2 + c1x/2)^2 = (c1^2/4-c2)x^2-c3x-c4 两边再同时加上 (x^2+c1x/2)y + (y/2)^2 (x^2 + c1x/2)^2 + (x^2+c1/2)y + (y/2)^2 = (c1^2/4-c2)x^2-c3x-c4 + (x^2+c1x/2)y + (y/2)^2 再整理： ((x^2+c1x/2）+ y/2)^2 = (c1^2/4-c2+y)x^2 + (c1y/2-c3)x +(y/2)^2-c4.....（1）左边为一个完全平方式，右边为一个x的二次多项式(系数中含有y), 为了使右边也成为一个完全平方式，必须使其判别式为零，即： DAT = ^2 - 4 * * = 0 于是，得到一个关于y的一元三次方程。 -y^3 + c2y^2 + (4c4-c1c3)y+c3^2+c4c1^2-4c4c2=0 只要解这个一元三次方程得到一个解y0，代入（1）右边，并令右边=0解出配方常量x0; 便可得到两边都为完全平方的方程，原一元四次方程被降为一元二次方程。可轻松求解。 ((x^2+c1x/2）+ y0/2)^2 =(x-x0)^2 即： (x^2+c1x/2）+ y0/2 = x-x0 (x^2+c1x/2）+ y0/2 = -x+x0 整理为： x^2+（c1/2-1）x + y0/2 + x0 = 0 x^2 +(c1/2+1）x + y0/2 - x0 = 0 解这两个一元二次方程的解，即为一元四次方程的解。原码核心部分如下： Procedure T1P4Func.Calculate; //一元四次方程求解计算方法 var Fa,Fb,Fc,Fd,Fe,C4 : Real; //系数计算中间变量 Za,Zb,Zc,zd,Ze : Real; //系数计算中间变量 AFunc : T1P3Func; //一元三次方程实例，用来计算配方系数 BFunc,CFunc,DFunc : T1P2Func; //一元二次方程实例，用来求解配方解和最终解 begin Fa := Coefficients ; //按通式ax^4+bx^3+cx^2+dx+e=0获取系数 Fb := Coefficients ; Fc := Coefficients ; Fd := Coefficients ; Fe := Coefficients ; Zb := Fb/Fa; //化简为通式：x^4+zbx^3+zcx^2+zdx+ze=0 Zc := Fc/Fa; Zd := Fd/Fa; Ze := Fe/Fa; Za := Zb*Zb; AFunc := T1P3Func.Create(3); //建立配方系数求解方程,通式：c1y^3+c2y^2+c3^y+c4=0; AFunc.Coefficients := -1; //-y^3 + Zc * y^2 + (4Ze-Zb*Ze)y + Zd*Zd + Ze*Zb*Zb - 4*Zc*Ze = 0 AFunc.Coefficients := Zc; AFunc.Coefficients := -Zb*Zd + 4*Ze; AFunc.Coefficients := Zd*Zd + Za*Ze-4*Zc*Ze; AFunc.Calculate; //求解配方系数y0 if AFunc.RootCount 0 then begin //将解出的y代入(c1^2/4-c2+y)x^2 + (c1y/2-c3)x +(y/2)^2-c4 = 0，以配出完全平方式 C4 := AFunc.RSolutions /2; //y0/2?????????这里可能需要判断，要取实根。 BFunc := T1P2Func.Create(2); //用配方解配方 BFunc.Coefficients := Za/4-Zc + AFunc.RSolutions ; BFunc.Coefficients := Zb*C4 - Zd; BFunc.Coefficients := C4*C4 - Ze; BFunc.Calculate; //求出配方的常数项x=x0,右边可配成(x-x0)^2 CFunc := T1P2Func.Create(2); //消去配方，化为两个对偶的一元二次方程 DFunc := T1P2Func.Create(2); // ((x^2+c1x/2）+ y/2)^2 =(x-x0)^2 CFunc.Coefficients := 1; //x^2+（c1/2-1）x + y0/2 + x0 = 0 CFunc.Coefficients := Zb/2-1; CFunc.Coefficients := C4+BFunc.RSolutions ; CFunc.Calculate; DFunc.Coefficients := 1; //x^2 +(c1/2+1）x + y0/2 - x0 = 0 DFunc.Coefficients := Zb/2+1; DFunc.Coefficients := C4-BFunc.RSolutions ; DFunc.Calculate; RSolutions := CFunc.RSolutions ;//解实部 RSolutions := CFunc.RSolutions ; RSolutions := DFunc.RSolutions ; RSolutions := DFunc.RSolutions ; VSolutions := CFunc.VSolutions ;//解虚部(如果有) VSolutions := CFunc.VSolutions ; VSolutions := DFunc.VSolutions ; VSolutions := DFunc.VSolutions ; FRootCount := 4; BFunc.Free; CFunc.Free; DFunc.Free; end else begin FRootCount := 0; end; AFunc.Free; end; 跟踪程序运行,常发现一元三次方程AFunc无解返回,导致求解失败. 进而查看一元三次方程求解程序问题. 一元三次方程求解用的是"盛金公式: 盛金公式　　一元三次方程aX^3＋bX^2＋cX＋d=0，（a，b，c，d∈R，且a≠0）。重根判别式： A=b^2－3ac； B=bc－9ad； C=c^2－3bd，总判别式： Δ=B^2－4AC。当A=B=0时，盛金公式① ： X1=X2=X3=－b/(3a)=－c/b=－3d/c。当Δ=B^2－4AC0时，盛金公式②： X1=(－b－((Y1)^(1/3)＋(Y2)^(1/3)))/(3a)； X2，3=(－2b＋(Y1)^(1/3)＋(Y2)^(1/3)±3^(1/2)((Y1)^(1/3)－(Y2)^(1/3))i)/(6a)，其中Y1，2=Ab＋3a(－B±(B^2－4AC)^(1/2))/2，i^2=－1。当Δ=B^2－4AC=0时，盛金公式③： X1=－b/a＋K； X2=X3=－K/2，其中K=B/A，(A≠0)。当Δ=B^2－4AC0时，盛金公式④： X1=(－b－2A^(1/2)cos(θ/3))/(3a)； X2，3=(－b＋A^(1/2)(cos(θ/3)±3^(1/2)sin(θ/3)))/(3a)，其中θ=arccosT，T= (2Ab－3aB)/(2A^(3/2))，(A0，－1T1)。关键原码如下: {T1P3Func} Procedure T1P3Func.Calculate; var Fa,Fb,Fc,Fd : Real; DA,DB,DC,DD : Real; Y1,Y2,K : Real; begin Fa := Coefficients ; Fb := Coefficients ; Fc := Coefficients ; Fd := Coefficients ; DA := Fb*Fb - 3*Fa*Fc; DB := Fb*Fc - 9*Fa*Fd; DC := Fc*Fc - 3*Fb*Fd; DD := DB*DB - 4 * DA * DC; if (DA = DB) and (DA =0) then begin RSolutions := -Fb/(3*Fa); RSolutions := -Fb/(3*Fa); RSolutions := -Fb/(3*Fa); VSolutions := 0; VSolutions := 0; VSolutions := 0; FRootCount := 1; end else if DD 0 then begin Y1:= DA * Fb + 3*Fa*(-Fb+sqrt(DD))/2; Y2:= DA * Fb + 3*Fa*(-Fb-sqrt(DD))/2; if (Y10) and (Y20) then begin Y1:= power(Y1,1/3); Y2:= power(Y2,1/3); end else if (Y10) and (Y20) then begin Y1:= power(Y1,1/3); Y2:= -power(-Y2,1/3); end else if (Y10) and (Y20) then begin Y1:= -power(-Y1,1/3); Y2:= -power(-Y2,1/3); end else if (Y10) and (Y20) then begin Y1:= -power(-Y1,1/3); Y2:= power(Y2,1/3); end; RSolutions := (-Fb - Y1 - Y2)/(3*Fa); VSolutions := 0; RSolutions := (-2*Fb+Y1+Y2)/(6*Fa); VSolutions := Sqrt(3)*(Y1-Y2)/(6*Fa); RSolutions := RSolutions ; VSolutions := -VSolutions ; FRootCount := 3; end else if DD = 0 then begin K := DB/DA; RSolutions := -Fb/Fa + K; RSolutions := -K/2; RSolutions := RSolutions ; VSolutions := 0; VSolutions := 0; VSolutions := 0; FRootCount := 3; end else if DD0 then begin K := (2*DA*Fb-3*Fa*DB)/(2*Power(DA,3/2)); Y1 := Arccos(K); RSolutions := (-Fb-2*Sqrt(DA)*Cos(Y1/3))/(3*Fa); RSolutions := (-Fb+Sqrt(DA)*(Cos(Y1/3)+ sqrt(3)*sin(Y1/3)))/(3*Fa); RSolutions := (-Fb+Sqrt(DA)*(Cos(Y1/3)- sqrt(3)*sin(Y1/3)))/(3*Fa); VSolutions := 0; VSolutions := 0; VSolutions := 0; FRootCount := 3; end; end; 仔细对照算法说明,并未发现程序实现和算法说明的差异.用实例测试发现: 这段程序在解DD0的情况时,不能得到正确的解. 如方程(x^2+1)(x+3)=0的解,一看就知道是:i,-i和-3. 展开为:x^3+3x^2+x+3=0 用程序解时:DD=12000(和手工验算符合) 得到的解是: -1.21822935060556+0i -0.890885324697222+2.13785617558042i -0.890885324697222-2.13785617558042i 和i,-i,-3相距甚远,难道盛金公式用的有问题?; 个人分类: 信息探索|13622 次阅读|4 个评论

更多...

帐号		自动登录	找回密码
密码			注册

关闭 安全验证

标签: 一元四次方程

相关帖子

相关日志

关闭安全验证