数学-林士谔算法


代数基本定理

1 代数基本定理

任何复系数一元n次多项式(n至少为1)方程在复数域上至少有一根。

n次复系数多项式方程在复数域内有且只有n个根,重根按重数计算。

证明不会

2 虚根成对定理

在实系数多项式分解中,虚根成对分解,实根单一分解,因此对于奇数次多项式,一定有实根。

简单理解: 假设一个根为虚根 /((x-a-bi)/) ,那么必须有另一个共轭根 /((x-a+bi)/) 这样,两者相乘可以得到实系数二次多项式 /(x^2-2ax+a^2+b^2/) ,否则无法消去虚数。

而且值得注意的是任何无实根的实系数多项式 /(x^2+px+q/) 都可以分解为/(x^2-2ax+a^2+b^2/) . (解出a,b从/(q>/frac{p^2}4/)入手证明)

3 林士谔算法

目标:找到多项式的根(包括实根和虚根)。

由代数基本定理,大于二次的多项式一定存在二次三项式因式。如何在实系数多项式中找到一个根,可以考虑在实系数多项式中找到一个二次三项式的因式。

因此可以用逐次分解二次三项式的方法找到所有根。

原理

/[f(x)=(x^2+px+q)g(x)
/]

我们要找到一个二次三项式,直接来找是很难的,可以先试着除一个二次三项式

/[f(x) = (x^2+px+q)g(x)+rx+s
/]

我们想办法通过修正p和q来消去余项。

/[p_1=p+dq
/]

/[q_1=q+dq
/]

r和s由p和q决定,对r和s求偏导,得到关系如下。

/[dr=/frac{/partial r}{/partial p}dp + /frac{/partial r}{/partial q}dq
/]

/[ds=/frac{/partial s}{/partial p}dp+/frac{/partial s}{/partial q}dq
/]

将f(x)分别对p、q求偏导,得到如下:

/[0=xg(x)+/frac{/partial g(x)}{/partial p}(x^2+px+q)+/frac{/partial r}{/partial p}x + /frac{/partial s}{/partial p}
/]

/[0=g(x)+/frac{/partial g(x)}{/partial q}(x^2+px+q)+/frac{/partial r}{/partial q}x + /frac{/partial s}{/partial q}
/]

由上二个式子可以变形为

/[xg(x)=-/frac{/partial g(x)}{/partial p}(x^2+px+q)-/frac{/partial r}{/partial p}x – /frac{/partial s}{/partial p}
/]

/[g(x)=-/frac{/partial g(x)}{/partial q}(x^2+px+q)-/frac{/partial r}{/partial q}x – /frac{/partial s}{/partial q}
/]

可以看到我们要求的偏导数正是xg(x)和g(x)对二次三项式的余式。可以直接计算。

因为我们希望r和s都修正到0,令dr=-r,ds=-s,可以近似方程为

/[-r=/frac{/partial r}{/partial p}dp + /frac{/partial r}{/partial q}dq
/]

/[-s=/frac{/partial s}{/partial p}dp+/frac{/partial s}{/partial q}dq
/]

可以近似解得dp和dq。

多次近似求解可以不断逼近。

代码如下

#include <bits/stdc++.h>
using namespace std;
const double eps = 1e-6;
double b[2333], c[2333];
void shie(double a[], int n, double &p, double &q)
{
    //b = g(x) = f(x) // (x^2 + px + q)
    memset(b, 0, sizeof(b));
    //c = x^2 * g(x) // (x^2 + px + q) from c[] we can get xg(x) % (x^2 + px + q) and g(x) % (x^2 + px + q)
    memset(c, 0, sizeof(c));
    p = 0, q = 0;
    double dp = 1;
    double dq = 1;
    while (dp > eps || dp < -eps || dq > eps || dq < -eps)
    {
        double p0 = p;
        double q0 = q;
        b[n - 2] = a[n];
        b[n - 3] = a[n - 1] - p0 * b[n - 2];
        c[n - 2] = b[n - 2];
        c[n - 3] = b[n - 3] - p0 * b[n - 2];
        for (int j = n - 4; j >= 0; --j)
        {
            b[j] = a[j + 2] - p0 * b[j + 1] - q0 * b[j + 2];
            c[j] = b[j] - p0 * c[j + 1] - q0 * c[j + 2];
        }
        double r = a[1] - p0 * b[0] - q0 * b[1];
        double s = a[0] - q0 * b[0];
        double rp = c[1];
        double sp = b[0] - q0 * c[2];
        double rq = c[0];
        double sq = -q0 * c[1];
        dp = (rp * s - r * sp) / (rp * sq - rq * sp);
        dq = (r * sq - rq * s) / (rp * sq - rq * sp);
        p += dp;
        q += dq;
    }
}
int main()
{
    double a[] = {-6, 11, -6, 1};
    int n = 3;
    double p, q;
    shie(a, n, p, q);
    cout << p << endl << q;
    return 0;
}

其中计算/(x^2g(x)/)对/((x^2+px+q)/)的除法可以得到所有想要结果,非常巧妙。

原创文章,作者:bd101bd101,如若转载,请注明出处:https://blog.ytso.com/277426.html

(0)
上一篇 2022年7月28日
下一篇 2022年7月28日

相关推荐

发表回复

登录后才能评论