数值计算讨论组

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


    [非]非线性方程加速方法几个算例

    分享

    Admin
    Admin

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

    [非]非线性方程加速方法几个算例

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

    实验四 非线性方程实根的加速法
    一、目的与要求:
    1. 掌握解非线性方程实根的简单迭代法的加速法、牛顿法加速法(下山法)与重根求解方法的上机编程运算.
    2. 通过上机编程运算比较不同方法的收敛速度.
    二、实验内容
    问题1. 用加权法加速技术求方程x=exp(-x)在0.5附近的一个根 .
    问题2. 用埃特金迭代法求方程x^2=x^3-1在x0=1.5附近的根.
    问题3. 用Steffenson算法求方程x^3 –x – 1=0在(1,1.5)内的根 .


    问题4. 分别用牛顿法与牛顿加速法(下山法)法求方程x^3 –x – 1=0在1.5附近的一个根.
    问题5. 已知x=sqrt(2)是方程f(x)=x^4-4x+4=0 的二重根,用牛顿切线法和重根修正公式求解.
    问题1
    代码:

    /*
    问题1 假设g(x)=exp(-x),加权加速方法为迭代 x(k+1) = g(x(k)), 修正 x(k+1) = (x(k+1)-L*x(k))/(1-L),其中L=dg/dx|(x_star),x_star = 0.5

    迭代4次就可以使得精度达到1e-4?

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

    int main()
    {
       int k;
       double x,xk,L,G,err;
       double g(double x);
       double dg(double x);

       k = 0;
       xk = 0.5;
       L = dg(xk);
       G = 1 - L;
       err = 1;

       while(err > 0.0001)
       {
          x = g(xk); // 迭代
          x = (x - L*xk) / G; // 修正
          err = fabs(x - xk);
          xk = x;
          k++;
       }

       cout << "K-th iteration: " << k << endl;
       cout << "err= " << setprecision(12) << err << endl;
       cout << "x= " << setprecision(12) << x << endl;
       cout << "g(x)= " << setprecision(12) << g(x) << endl;
       return 0;
    }

    double g(double x)
    {
       return (exp(-x));
    }

    double dg(double x)
    {
       return (-exp(-x));
    }

    问题2
    代码:

    /*
    问题2 用Aitken Method迭代法求方程x^2=x^3-1在x0=1.5附近的根. (根所在的区间为[1,2])
    方程等价于 x = (x^2+1)^(1/3),简单迭代法适用其上
    Aitken Method:x(k+3) = x(k) - (x(k+1)-x(k))^2/(x(k+2)-2x(k+1)+x(k)),首两项可用简单迭代法给出
    */
    #include <iostream>
    #include <cmath>
    #include <iomanip>
    using namespace std;

    int main()
    {
       int k;
       double x,xk0,xk1,xk2,err;
       double f(double x);
       
       k = 0;
       xk0 = 1.5;
       xk1 = f(xk0);
       xk2 = f(xk1);
       err = 1.0;
       while(err > 0.00001)
       {
          // 迭代序列 xk0,xk1,xk2,x...
          x = xk0 - pow((xk1 - xk0),2)/(xk2 - 2*xk1 + xk0);
          xk0 = xk1;
          xk1 = xk2;
          xk2 = x;
          err = fabs(xk2 - xk1);
          k++;
       }

       cout << "K-th iteration: " << k << endl;
       cout << "err= " << setprecision(12) << err << endl;
       cout << "x= " << setprecision(12) << x << endl;
       cout << "f(x)= " << setprecision(12) << f(x) << endl;
       return 0;
    }

    double f(double x)
    {
       return cbrt(pow(x,2) + 1.0);  // cbrt 立方根函数
    }
    问题3
    代码:

    /*
    问题3 用Steffenson算法求方程x^3 –x – 1=0在(1,1.5)内的根
    Steffenson Method:yk = f(x(k)), zk = f(yk), x(k+1) = x(k) - (yk - x(k))^2/(zk - 2yk + x(k));
    */
    #include <iostream>
    #include <cmath>
    #include <iomanip>
    using namespace std;

    int main()
    {
       int k;
       double x,xk,yk,zk,err;
       double f(double x);
       
       k = 0;
       xk = 1.5;
       yk = f(xk);
       zk = f(yk);
       err = 1.0;
       while(err > 0.0001)
       {
          // 迭代序列 xk,x...
          x = xk - pow((yk - xk),2)/(zk - 2*yk + xk);
          yk = f(x);
          zk = f(yk);
          err = fabs(x - xk);
          xk = x;
          k++;
       }

       cout << "K-th iteration: " << k << endl;
       cout << "err= " << setprecision(12) << err << endl;
       cout << "x= " << setprecision(12) << x << endl;
       cout << "f(x)= " << setprecision(12) << f(x) << endl;
       return 0;
    }

    double f(double x)
    {
       return (pow(x,3) - 1.0); 
    }
    问题4
    代码:

    /*
    问题4 分别用牛顿法与牛顿加速法(下山法)法求方程x^3 –x – 1=0在1.5附近的一个根
    f(x) = x^3 - x - 1
    仅考虑牛顿加速法(下山法):x(k+1) = x(k) - tk * f(x(k))/f '(x(k))
    初值选为0.6
    */

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

    int main()
    {
       int k;
       double x,xk,fxk,t,direction,err;
       double f(double x);
       double df(double x);
       
       k = 0;
       xk = 0.6;
       err = 1.0;
       while(err > 0.000001)
       {
          // 迭代序列 xk,x...
          t = 1.0;
          fxk = f(xk);
          direction = fxk/df(xk);
          while ( true )
          {
             x = xk - t*direction;
             if (fabs(f(x)) < fabs(fxk))
                break;
             else
             {
                t = t/2.0;
                continue;
             }
          }      
          err = fabs(x - xk);
          xk = x;
          k++;
          // 每步迭代过程
          cout << "K-th iteration: " << k << endl;
          cout << "t = " << t << endl;
          cout << "err= " << setprecision(12) << err << endl;
          cout << "x= " << setprecision(12) << x << endl;
          cout << "f(x)= " << setprecision(12) << f(x) << endl;
          cout << endl;
       }

       return 0;
    }

    double f(double x)
    {
       return (pow(x,3) - x - 1.0);
    }
    double df(double x)
    {
       return (3.0*pow(x,2) - 1.0);
    }
    问题5
    代码:

    /*
    问题5. 已知x=sqrt(2)是方程f(x)=x^4-4*x^2+4=0 的二重根,用牛顿切线法和重根修正公式求解
    仅考虑重根修正公式:x(k+1) = x(k) - m*f(x(k))/f ' (x(k)), m为根的重数
    */
    #include <iostream>
    #include <cmath>
    #include <iomanip>
    using namespace std;

    int main()
    {
       int k;
       double x,xk,m,err;
       double f(double x);
       double df(double x);
       
       k = 0;
       m = 2.0;
       xk = 1.5;
       err = 1.0;
       while (err > 0.0001)
       {
          // 迭代序列 xk,x...
          x = xk - m*f(xk)/df(xk);
          err = fabs(x - xk);
          xk = x;
          k++;
       }

       cout << "K-th iteration: " << k << endl;
       cout << "err= " << setprecision(12) << err << endl;
       cout << "x= " << setprecision(12) << x << endl;
       cout << "f(x)= " << setprecision(12) << f(x) << endl;
       return 0;
    }

    double f(double x)
    {
       return (pow(x,4) - 4.0*pow(x,2) + 4.0);
    }

    double df(double x)
    {
       return (4.0*pow(x,3) - 8.0*x);
    }

      目前的日期/时间是周五 十二月 14, 2018 2:32 am