yuange1975 发表于 2023-8-21 14:57:21

求任意多项式y=f(x)/g(x)的极值方程 带验证代码了

求任意多项式y=f(x)/g(x)的极值方程

经常会遇到求多项式y=f(x)/g(x)的极值,其中f和g是多项式。有时候这个函数的导数方程 y’=(f’*g-f*g’)/g^2=0 是高次方程又不能分解因数,这时候极值点要么没有解析表达式要么表达式非常复杂,再求原式的极值,就非常困难或者表达式非常复杂。

这时候我们就可以改变思路,去求极值满足的方程h(y)=0。

可以用程序来算,算法都比较简单。就是一些数组加减乘除和矩阵的计算。

y’=(f’*g-f*g’)/g^2

所以等式方程就是:

f’*g-f*g’=0,

等式方程次数n,一般情况次数n=f次数+g次数-1。

多项式加减乘除和幂方,就是数组的计算都很简单。计算后消去大于等于n的高次项。最后可以用一个n纬的数组表示,对应x的0到n-1次方的系数。

多项式的除法是先计算多项式的逆,就是已知g,求出满足g*G=1=I0的G。这样1/g=G。然后就变成乘法计算了。多项式的逆计算比加减乘除和幂方稍微复杂一点,需要计算一个矩阵的逆。

多项式消去高次项算法,就是类似消元法,从高次一步步减就行了,直到次数小于n。

多项式的逆,先消去高次。如果得到的次数大于0,就先求出g的逆多项式,就是求一个矩阵G的逆矩阵1/G,乘以多项式1记为I0。矩阵G的每行数组就是g乘以I0到In-1(就是对应x^(n-1))再消去高次项后的数组。

y=f*(1/G*I0)再消去高次项得到y=F,这时候求h(y)需要再求一个逆矩阵A的逆1/A。这个矩阵A由F的0次幂到n-1次幂消去高次项后到数组构成。

最后h(y)就是方程:

y^n-1/A*F^n=0


大致算法是这样,可能会有中间一些表诉的地方有小错误。等有空写出来程序代码,自动计算,只需要输入f和g的系数,实现快速程序计算输出方程 y^n-1/A*F^n=0 。写出代码后这中间的表诉小错误或者表诉不准确的地方也就能够清楚了。


程序出来后就不会得出图中这种复杂的极值表达式。比如图中y’=0是4次方程,极值方程也就是4次方程,极值表达式也和极值点一样的4次方程的根的形式,会比下面那个表达式简单简化很多。

并且得到极值方程后,再迭代得到数值解或者判断界限就是得到一些缩放的不等式等都很方便。

如果y'=0能分解,先直接分解计算简单一些,如果程序计算就无所谓了。y'=0能分解的话,不分解得到的h(y)=0也是一样的能分解。

这个也算是多项式的剩余环的计算了。




验证程序来了,使用Python的符号运算库整个过程就非常简单。验证代码暂时没有考虑任意函数任意次数的变量个数等细节问题。

运算结果得到方程:
y**4 - 45*y**3 + 2817*y**2/4 - 4860*y + 13068=0

其最大根21.2866399744823和后面的图像对应得起来,还有最小根9.58346869275525也是对应得起来的。



Python 代码如下:
# y=f(x)/g(x) 的极值方程h(y)=0 求解器
# copy by yuange 2023.8.20
from sympy import *
x,y=symbols('x y')
a,b,c,d=symbols('a b c d')
A,B,C,D=symbols('A B C D')
f=2*(2*x**2+3)*(3*x**2+2*x+3)
g=(x*x+1)**2
gg=a*x**3+b*x**2+c*x+d
fg=f/g
fg1=cancel(diff(fg))
fg2=fraction(fg1)
yy=fg2
eq0=Eq(1,div(g*gg,yy))
res0=solve(eq0,)
g1=gg.subs(res0)
print(g1)
ff=div(f*g1,yy)
hy=ff**4+A*ff**3+B*ff*ff+C*ff+D
hy2=y**4+A*y**3+B*y**2+C*y+D
eq1=Eq(0,div(hy,yy))
res=solve(eq1,)
hy3=hy2.subs(res)
print(hy3)
res2=nsolve(hy3,20)
print(res2)

页: [1]
查看完整版本: 求任意多项式y=f(x)/g(x)的极值方程 带验证代码了