数值计算讨论组

欢迎来到数值计算讨论组。 本论坛旨在增强师生交流的时效性,主要讨论数值计算相关问题,包括初等数值分析、矩阵计算、非线性方程(组)数值解法、微分方程数值解法等方面的内容。


    [非]多项式求根劈因子法

    分享

    Admin
    Admin

    帖子数 : 46
    注册日期 : 10-05-13

    [非]多项式求根劈因子法

    帖子 由 Admin 于 周六 五月 15, 2010 4:47 pm

    该问题不是实验大纲中的问题。

    问题:用劈因子法求f(x) = x^6 + 8*x^5 + 40*x^4 + 70*x^3 + 89*x^2 + 62*x + 50 = 0的一个二次因式,该二次因式有x^2 + u*x + v的形式,进而由x^2 + u*x + v可以得到原多项式方程的两个共轭复根。程序中一些符号的含义与课件中一致。

    代码:

    #include <iostream>
    #include <cmath>
    #include <iomanip>
    using namespace std;

    int main()
    {
       int k, i;
       const int n = 6;
       double u,v,err;
       double r0,r1,s0,s1;
       double D,D1,D2,delta_u,delta_v;
       double a[n+1] = {1.0, 8.0, 40.0, 70.0, 89.0, 62.0, 50.0};
       double b[n-1],c[n-1];

       k = 0;
       err = 1.0;
       u = a[n-1]/a[n-2];  // u,v的初值
       v = a[n]/a[n-2];
       b[0] = a[0];
       c[0] = b[0];
       while(err > 0.000001)
       {
          b[1] = a[1] - u*b[0];
          for (i = 2; i<=n-2; i++)
          {
             b[i] = a[i] - u*b[i-1] - v*b[i-2];
          }
          r0 = a[n-1] - u*b[n-2] - v*b[n-3];
          r1 = a[n] - v*b[n-2];

          c[1] = b[1] - u*b[0];
          for (i = 2; i<=n-2; i++)
          {
             c[i] = b[i] - u*c[i-1] - v*c[i-2];
          }
          s0 = r0 - u*c[n-2] - v*c[n-3];
          s1 = c[n-2] + u*c[n-3];
    /*
          每次循环都要求解一个2元1次方程组,未知量分别为u,v的修正delta_u,delta_v         
           /  r0 + dr0/du * x + dr1/du * y = 0
          |
           \  r1 + dr1/dv * x + dr1/dv * y = 0         
          其中delta_u = x, delta_v = y
          dr0/dv = -s0,        dr1/dv = -s1
          dr0/du = u*s0 - s1,  dr1/du = v*s0
          求解该方程组直接用了Cramer法则
    */
          D = -(u*s0 - s1)*s1 + v*s0*s0;
          D1 = r0*s1 - r1*s0;
          D2 = -(u*s0 - s1)*r1 + v*s0*r0;
          delta_u = D1/D;
          delta_v = D2/D;

          u = u + delta_u;
          v = v + delta_v;
           err = sqrt(pow(delta_u,2)+pow(delta_v,2));
          k++;
       }

       cout << "K-th iteration: " << k << endl;
       cout << "err= " << setprecision(12) << err << endl;
       cout << "u= " << setprecision(12) << u << endl;
       cout << "v= " << setprecision(12) << v << endl;
       return 0;
    }
    // n次多项式f(x)
    double f(double x)
    {
       // 计算x的幂可以比下面的方式更简单
       return (pow(x,6) + 8.0*pow(x,5) + 40.0*pow(x,4) + 70.0*pow(x,3) + 89.0*pow(x,2) + 62.0*x + 50.0);
    }
     

      目前的日期/时间是周一 五月 21, 2018 5:04 pm