数学建模社区-数学中国

标题: 用 Mathematica 求解多项式 [打印本页]

作者: Asdmath2020 时间: 2020-3-17 14:54
标题: 用 Mathematica 求解多项式

多项式是由一组常数系数，a、b、c、……（数值）确定的。
TableForm[{a x + b, a x^2 + b x + c, a x^3 + b x^2 + c x + d, ". . ."}] //

多项式求解问题就是找到一个值 x，使这些项的总和等于 0. 根据 x 的最高次数分别称为线性、二次、三次、四次、五次、六次、七次、八次...... 多项式。我们称 y = a x + b 为线性，是因为它的图线是一条直线. 比如令 a = 2，b = 3，

Plot[2 x + 3, {x, -2, 1}]

2 x + 3 = 0 的（唯一）解是 x = -3/2. 一般情况下，有 x = -b/a. 由于含有 x 的平方项，y = a x^2 + b x + c 是二次的. 你会记得一元二次方程有两个通解：

Solve[a x^2 + b x + c == 0, x]

这样的表达式被称为不尽根式. 最常见的应用是在几何上. 圆、抛物线和双曲线通常由二次多项式指定。当我们想知道一个二次多项式与已知直线何时相交时，我们就得到一个二次方程. 这甚至发生在双曲线是矩形的，例如：
Plot[{(1 - x)/(x + 2), 2 x + 3}, {x, -4, 1}, PlotRange -> {{-4, 1}, {-6, 6}},

表示双曲线的分支与直线相交的方程为

(1 - x)/(x + 2) == 2 x + 3;

要得到形如 x^2 + b x + c ==0 的方程，需要两边同时乘以 x + 2，

Distribute[(x + 2)*%, Equal][color=rgba(0, 0, 0, 0.498)]1 - x == (2 + x) (3 + 2 x)

再两边同时减去(2 + x) (3 + 2 x)，

(2 + x) (3 + 2 x) - # & /@ %

-1 + x + (2 + x) (3 + 2 x) == 0

然后展开：

Expand@%

5 + 8 x + 2 x^2 == 0
5 O1 _; \ G/ j# N7 L* Y

因为 Mathematica 不是使用传统的方法处理方程，"两边同时减去"这一步需要解释一下. 它使用 & 定义了一个函数，从 (2 + x) (3 + 2 x)减去（#），然后"映射" ( /@ ) 到等式两边.
不管怎样，求解二次方程：
Solve[%]
{{x -> 1/2 (-4 - Sqrt[6])}, {x -> 1/2 (-4 + Sqrt[6])}}
或者近似; U! U; B* H1 i3 R6 p( z0 N( H

N@%
{{x -> -3.22474}, {x -> -0.775255}}

- Y6 F3 ~( p0 y7 V) W' w
当然，Mathematica 非常乐意直接求解：
Solve[(1 - x)/(x + 2) == 2 x + 3]
{{x -> 1/2 (-4 - Sqrt[6])}, {x -> 1/2 (-4 + Sqrt[6])}}
我们可以把这些值作为绘图范围
' a/ j: I) p) g
Flatten[% /. Rule -> List]
{x, 1/2 (-4 - Sqrt[6]), x, 1/2 (-4 + Sqrt[6])}
来直观地检查合理性：
4 {* r6 c& P! j( F; l5 y
Plot[{(1 - x)/(x + 2), 2 x + 3}, %[[{1, 2, 4}]]]# e' ~# R* i% m6 a* {- E

2020-3-17 14:32 上传
下载附件 (8.56 KB)

该图形以交点作为起点和终点. （ {1, 2, 4} 是为了去除额外的 x.）

要是双曲线是颠倒的，貌似就没解了：

Plot[{(x - 1)/(x + 2), 2 x + 3}, {x, -4, 1}, PlotRange -> {{-4, 1}, {-6, 6}},
AspectRatio -> 4/5]

但 Mathematica 很聪明，它给的结果不是{}，而是

Solve[(x - 1)/(x + 2) == 2 x + 3]

{{x -> 1/2 (-3 - I Sqrt[5])}, {x -> 1/2 (-3 + I Sqrt[5])}}

也就是说，解含有虚数. 二次方程是通过配方法来求解的，两边同时加上b^2/4a-c:

# + b^2/(4 a) - c & /@ (a x^2 + b x + c == 0)

b^2/(4 a) + b x + a x^2 == b^2/(4 a) - c

然后左边配成平方除以4a：

Factor /@ %

(b + 2 a x)^2/(4 a) == -((-b^2 + 4 a c)/(4 a))

然后乘以 4a 后

4 a # & /@ %

(b + 2 a x)^2 == b^2 - 4 a c

我们可以取两边的平方根.

Sqrt /@ % // PowerExpand

b + 2 a x == Sqrt[b^2 - 4 a c]

现在它变成了一个线性方程. PowerExpand 强制 Mathematica 选择 b + 2 a x，而不是 - b - 2 a x 作为平方的方根.

正如"每个人"都记着二次方程的解，"没有人"记得三次方程的解。原因是：

Solve[ a x^3 + b x^2 + c x + d == 0, x] // TraditionalForm

现在考虑一个简单情形：
Plot[x^3 - x + 1/3, {x, -3/2, 3/2}]

显然它有三个实根，在 - 1.1, .35, 和 .75 附近. 具体是多少呢？

Solve[x^3 - x + 1/3 == 0]

嘿，等一下，根不是实数吗？是不是 Mathematica 忘了化简了？

Simplify@%%

这么复杂的式子怎么会是实数呢？让我们看一下虚部：

ComplexExpand@Im[x /. %]

{0, 0, 0}

那好吧，老天爷，告诉我们实部是啥吧！

ComplexExpand@Re[x /. %%]

三角函数？！还不如平方根和立方根呢！

Developer`TrigToRadicals@%%

天啊！那些三角函数是实数，但为啥这里却跟着一堆虚数单位？- 1.1、.35 和 .75 在哪呢？？请给出数来.

N@%

{0.742227 + 0. I, 0.394931 + 0. I, -1.13716 + 0. I}

它们就在那，但这些 0. i 是怎么回事？奇怪的是，它们是不可避免的. 作为一个数学分支的伽罗瓦理论已经证明，不含虚数立方根的解式是不存在的，即使它们的加和为实数.

Henry Baker 的动画（本文顶部可以看到实际动画）展示了均为实数的三个根之间的关系：

这是张一般情况下的图片——三个实根的三次方程有一个拐点，它们关于拐点对称. 如果将拐点平移到原点，则会得到一个奇函数 f（-x ）= -f（x）.

四次方程可以通过将两条曲线相交得到. 一般情况下的四次方程会让人有点抓狂，如果不怕的话就按住 shift return 键试试吧.

Solve[ a x^4 + b x^3 + c x^2 + d x + e == 0, x]

求解五次方程就更是不要命了.

Solve[ a x^5 + b x^4 + c x^3 + d x^2 + e x + f == 0, x]

哈哈，Mathematica 放了我们一马，结果被剪切了，但为什么不至少提供一下互动大型表达式浏览器呢？因为它不能. 不存在一般五次方程的根表达式. 显然，通过因式分解我们可以求解某些五次方程.

Expand@# == # &[(x^2 + 1) (x^3 - x - 1)]

-1 - x - x^2 + x^5 == (1 + x^2) (-1 - x + x^3)

历经几个世纪的挫败，求解五次方程已经与三等分角和倍立方问题一样成为困扰人们的几大数学难题。

许多人错误地认为唯一可解的五次方程要么是可因式分解的，要么是显而易见的，如（x + a）^ 5 + b = 0. 但只有一小部分，接近0％，可以巧妙地解决，比如：

x^5 - 5 x^2 - 3 == 0

它的根相当繁琐. 唯一的实根是

至少从外观看来是实数. 平凡的 Mathematica 无法求解——连验证都不行！只能近似：

%^5 - 5`69 %^2

3.0000000000000000000000000000000000000000000000000000000000000000000

我曾经（正确地？）说服自己，每个可求解的六次方程都可以降次到具有二次不尽根系数的三次方程或具有三次不尽根系数的二次方程. 但谁会想要求解这样一个方程呢？几何再次派上用场了. 问题：将一个正方形拆分成有限个锐角、等腰三角形. 可以用十个：

确定 Subscript[x, 1]、Subscript[y, 1]、Subscript[x, 2] 和 Subscript[r, 3] 的方程是12次的！这些方程可以通过因式分解降为具有不尽根系数的六次方程，但专家 Noam Elkies 认为这是不可求解的. 然而当次数为6、8、9、……或任何复合数（或称非质数）时，有时候会比较幸运.

Factor[5 + 4 x^2 - 4 x^3 + x^4 - 2 x^5 + x^6]

5 + 4 x^2 - 4 x^3 + x^4 - 2 x^5 + x^6

不能进行因式分解. 然而

Solve[% == 0]

可以求得所有六个解！怎么实现呢？这个六次方程可以写作

x^2 - 4 x + 5 /. x -> x^3 - x^2

5 - 4 (-x^2 + x^3) + (-x^2 + x^3)^2

Factor@%

5 + 4 x^2 - 4 x^3 + x^4 - 2 x^5 + x^6

也就是说，将三次方程替换为二次方程。如果我们注意到这一点，我们只是用y来代替 x ^ 3 - x ^ 2 - 2 ，对得到的二次方程求 y，然后求解关于 x 的三次方程，用 y 表示。我们是怎么注意到这一点的？用魔法函数

Decompose[5 + 4 x^2 - 4 x^3 + x^4 - 2 x^5 + x^6, x]

{5 - 4 x + x^2, -x^2 + x^3}

这正好是 Solve 函数所了解的. 知道吗？你的八次方程可能只是三个二次方程的组合.

但请注意：这个六次方程的解，既不能因式分解

5 + 18 x + 36 x^2 + 36 x^6

也不能正常分解.

{570630428688384000000 + 4891824455002619904 x + 161093791317491712
x^2 + 12153384861696 x^3 + 984379392 x^4 - 17280 x^5 + x^6}

Simplify[% /. x -> %%]

{0}

所以括号中的数量满足最小次数为36的最小多项式！

令人惊奇的是，这是一个甚至连 Mathematica 第11版都不知道的诀窍：如果系数形成回文，六次甚至八次方程总是可以求解！例如，

Solve[1 + x + x x - x^3 + x^4 + x^5 + x^6 == 0]

（失败。）但是我们可以对付完全一般的情况！对于任意a，b，c，d，假设 x 满足互反多项式（所以被称作回文式）

a x^6 + b x^5 + c x^4 + d x^3 + c x x + b x + a

a + b x + c x^2 + d x^3 + c x^4 + b x^5 + a x^6

现在假设 y = x + 1/x （或写作 x y = x^2 + 1），求六次多项式除以这个二次多项式（关于 x）的余式：

Factor[PolynomialRemainder[%, x + 1/x - y, x]]

(-x - y + x y^2) (-2 b + d - 3 a y + c y + b y^2 + a y^3)

求这个余式意味着减去二次式的倍数，使得六次多项式将简化为关于 x 的线性多项式. 但是我们假定了二次和六次多项式都是0，所以我们从0减去0，得到x和y之间的可疑关系，乘以我们可以求解的 y 的三次式！通过 y = x + 1 / x 来求解x.

回文多项式被称作互逆多项式的原因是，如果用 1/x 代替 x，两者具有相同的根，从而将系数的次序逆转（并除以 x^6）.

这个令 y = x + 1/x 的技巧可以成功的关键是我们可以将各项用它们的倒数匹配，并利用关系：

1/x^3 + x^3 == -3 (1/x + x) + (1/x + x)^3

1/x^2 + x^2 == -2 + (1/x + x)^2

朱利安和我有一个七次方程求解程序，但他不相信它能找到所有的解. 超过七次以后，能找到一个强有力的求解器机会会大大减小，TA在理论上可以求解的概率也是如此. 但是如果你的问题不是随机组成的，那么总是值得一试.

若想咨询Mathematica软件，欢迎联系武汉墨光科技market@asdoptics.com

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)