计算科学 - 使用 BDF 和 RK4 方法在 C++ 中求解此 ODE 耦合系统 - 吾爱随笔录

使用 BDF 和 RK4 方法在 C++ 中求解此 ODE 耦合系统

计算科学颂数字 C++ 稳定龙格库塔

2021-12-24 02:40:13

我正在尝试使用 BDF 4 阶方法解决 ODE 系统。我使用 RK4 找到前 3 个点，然后对于 BDF 的隐式部分，我使用 Newton-Raphson 迭代。不幸的是，我的解决方案失败了，但我不知道我做错了什么。我认为这与牛顿迭代有关，因为我对那部分不太自信。ODE 的形式为

$y_1' = 40·y_1 -2·y_2^2 +40·y_2^2 -100·y_1·y_2^2+...$ 大约还有 30 个术语

$y_2' = 2·y_1 +40·y_2 +...$ 大约还有 30 个术语

到目前为止，这是我的代码

#include<stdio.h>
#include<math.h>
#include<iostream>

#define MAX_N 10000
using namespace std;
static   double F(double, double);
static   double S(double, double);
static   double Fprime(double, double);
static   double Sprime(double, double);

int main()
{
    double K1_1,K2_1,K3_1,K4_1,K1_2,K2_2,K3_2,K4_2,W[MAX_N],V[MAX_N],H,T;
    double A = 0.0;
    double B = 5.0;
    int N = 10000;
    int I;

    cout.setf(ios::fixed,ios::floatfield);
    cout.precision(9);

    H = (B - A) / N;
    T = A;
    W[0] = 0.5;
    V[0] = 0.0;

    for (I=1; I<=3; I++)
    {
         K1_1 = H*F(W[I-1], V[I-1]);
         K1_2 = H*S(W[I-1], V[I-1]);

         K2_1 = H*F(W[I-1] + K1_1/2.0, V[I-1] + K1_2/2.0);
         K2_2 = H*S(W[I-1] + K1_1/2.0, V[I-1] + K1_2/2.0);

         K3_1 = H*F(W[I-1] + K2_1/2.0, V[I-1] + K2_2/2.0);
         K3_2 = H*S(W[I-1] + K2_1/2.0, V[I-1] + K2_2/2.0);

         K4_1 = H*F(W[I-1] + K3_1, V[I-1] + K3_2);
         K4_2 = H*S(W[I-1] + K3_1, V[I-1] + K3_2);

         W[I] = W[I-1] + 1/6.0*(K1_1 + 2.0*K2_1 + 2.0*K3_1 + K4_1);
         V[I] = V[I-1] + 1/6.0*(K1_2 + 2.0*K2_2 + 2.0*K3_2 + K4_2);

         T = A + I * H;

         cout <<"At time "<< T <<" the solution = "<< W[I] << endl;
    }
    //BDF order 4 to get the rest of the points
      for(I = 4; I <= N; I++)
      {
          //Newton Raphson method to get the values of W[I],V[I] for the implicit BDF
          double W_temp = W[I-1];
          double V_temp = V[I-1];
          double tol = 1e-14;
          double error = tol + 1;
          int iteration = 0;

          //Checking tolerance, the denominator not being too small, and a reasonable number of iterations
          while (error > tol && fabs(Fprime(W_temp, V[I-1]))>1e-14 && iteration < 1000)
          {
            W[I] = W_temp - F(W_temp, V[I-1])/Fprime(W_temp, V[I-1]);

            error = fabs(W[I] - W_temp);
            W_temp = W[I];
            iteration++;

          }
          iteration = 0;

          while (error > tol && Sprime(W[I-1], V_temp)>1e-14 && iteration < 1000)
          {
            V[I] = V_temp - S(W[I-1], V_temp)/Sprime(W[I-1], V_temp);

            error = fabs(V[I] - V_temp);
            V_temp = V[I];
            iteration++;

          }

          //BDF order 4
          W[I] = (48.0*W[I-1] - 36.0*W[I-2] + 16.0*W[I-3] - 3.0*W[I-4] + 12.0*H*F(W[I],V[I]))/25.0;
          V[I] = (48.0*V[I-1] - 36.0*V[I-2] + 16.0*V[I-3] - 3.0*V[I-4] + 12.0*H*S(W[I],V[I]))/25.0;

          T = A + I * H;

          cout <<"At time "<< T <<" the solution = "<< W[I] << endl;
      }

    return 0;
}

/*  First incremental function  */
double F(double y1, double y2)
{
   double f; 

   f = 40.0*y1 - 2.0*y2 + 40.0*pow(y2,2) - 100.0*y1*pow(y2,2) + 160.0*pow(y1,2)*pow(y2,4) + 100.0*pow(y1,2)*pow(y2,2) - 180.0*y1*pow(y2,6) + 180.0*y1*pow(y2,4) -240.0*pow(y1,4)*pow(y2,2) +100.0*pow(y1,3)*pow(y2,4) + 220.0*pow(y1,2)*pow(y2,6) - 180.0*y1*pow(y2,8) + 4.0*y1*y2 - 60*pow(y2,12) - 20.0*pow(y1,7) - 20.0*pow(y2,14) - 60.0*pow(y1,5)*pow(y2,2) + 180.0*pow(y1,4)*pow(y2,4) - 120.0*pow(y1,3)*pow(y2,6) - 120.0*pow(y1,2)*pow(y2,8) + 180.0*y1*pow(y2,10) + 100.0*pow(y1,6)*pow(y2,2) - 180.0*pow(y1,5)*pow(y2,4) + 100.0*pow(y1,4)*pow(y2,6) + 100.0*pow(y1,3)*pow(y2,8) - 180.0*pow(y1,2)*pow(y2,10) + 100.0*y1*pow(y2,12) + 140.0*pow(y2,8) - 20.0*pow(y2,6) - 100.0*pow(y1,3) + 80.0*pow(y1,5) + 20.0*pow(y2,10) - 4.0*pow(y2,3) - 100.0*pow(y2,4);
   return f;
}
double Fprime(double y1, double y2)
{
  double fprime;

  fprime = 40.0 - 100.0*pow(y2,2) + 480.0*pow(y1,2)*pow(y2,2) - 320.0*y1*pow(y2,4) + 200.0*y1*pow(y2,2) - 180.0*pow(y2,6) + 180.0*pow(y2,4) - 960.0*pow(y1,3)*pow(y2,2) + 300.0*pow(y1,2)*pow(y2,4) + 440.0*y1*pow(y2,6) - 180.0*pow(y2,8) + 4.0*y2 - 140.0*pow(y1,6) - 300.0*pow(y1,4)*pow(y2,2) + 720.0*pow(y1,3)*pow(y2,4) - 360.0*pow(y1,2)*pow(y2,6) - 240.0*y1*pow(y2,8) + 180.0*pow(y2,10) + 600.0*pow(y1,5)*pow(y2,2) - 900.0*pow(y1,4)*pow(y2,4) + 400.0*pow(y1,3)*pow(y2,6) + 300.0*pow(y1,2)*pow(y2,8) - 360.0*y1*pow(y2,10) + 100.0*pow(y2,12) - 300.0*pow(y1,2) + 400.0*pow(y1,4);
  return fprime;
}



/*  Second incremental function  */
double S(double y1, double y2)
{
   double s; 

   s = 2.0*y1 + 40.0*y2 - 2.0*pow(y2,2) - 100.0*pow(y1,2)*y2 + 200.0*y1*pow(y2,3) - 320.0*pow(y1,3)*pow(y2,3) + 420.0*pow(y1,2)*pow(y2,5) + 160.0*pow(y1,2)*pow(y2,3) - 200.0*y1*pow(y2,7) - 320.0*y1*pow(y2,5) - 60.0*pow(y1,4)*pow(y2,3) + 240.0*pow(y1,3)*pow(y2,5) - 360.0*pow(y1,2)*pow(y2,7) + 240.0*y1*pow(y2,9) + 120.0*pow(y1,5)*pow(y2,3) - 300.0*pow(y1,4)*pow(y2,5) + 400.0*pow(y1,3)*pow(y2,7) - 300.0*pow(y1,2)*pow(y2,9) + 120.0*y1*pow(y2,11) - 20.0*pow(y1,6)*y2 + 80.0*pow(y1,4)*y2 - 100.0*pow(y2,3) - 20.0*pow(y2,13) + 20.0*pow(y2,9) + 140.0*pow(y2,7) - 60.0*pow(y2,11) - 20.0*pow(y2,5);
   return s;
}
double Sprime(double y1, double y2)
{
  double fprime;

  fprime = 40.0 - 4.0*y2 - 100.0*pow(y1,2) + 600.0*y1*pow(y2,2) - 960.0*pow(y1,3)*pow(y2,2) + 2100.0*pow(y1,2)*pow(y2,4) + 480.0*pow(y1,2)*pow(y2,2) - 1400.0*y1*pow(y2,6) - 1600.0*y1*pow(y2,4) - 180*pow(y1,4)*pow(y2,2) + 1200.0*pow(y1,3)*pow(y2,4) - 2520.0*pow(y1,2)*pow(y2,6) + 2160.0*y1*pow(y2,8) + 360*pow(y1,5)*pow(y2,2) - 1500.0*pow(y1,4)*pow(y2,4) + 2800.0*pow(y1,3)*pow(y2,6) - 2700.0*pow(y1,2)*pow(y2,8) + 1320.0*y1*pow(y2,11) - 20.0*pow(y1,6) + 80.0*pow(y1,4) - 300.0*pow(y2,2) - 260.0*pow(y2,12) + 180.0*pow(y2,8) + 980.0*pow(y2,6) - 660.0*pow(y2,10) - 100.0*pow(y2,4);
  return fprime;
}

2个回答

从你的代码中，有几件事作为潜在的问题突然出现在我身上。

您似乎正在单独使用牛顿法 $y_1$ 和 $y_2$ 变量。这与对涉及两者的非线性系统使用牛顿法不同 $y_1$ 和 $y_2$ . 对于完整的牛顿法，您将在每一步求解一个 2x2 线性系统。你不仅要计算 $\partial F/\partial y_1$ 和 $\partial S/\partial y_2$ 但是也 $\partial F/\partial y_2$ 和 $\partial S/\partial y_1$ . 您写下的是一种坐标下降类型的方法，它可以解决一些问题，但这不是牛顿的方法。
牛顿方法不能保证收敛于你扔给它的任何初始数据。如果最初的猜测在二次收敛盆地之外很远，迭代可能会飞到无穷大。通常，您需要像线搜索或信任区域方法这样的全球化策略。
您的问题中的系数足够大，以至于在您解决的时间间隔内，如果您遇到溢出或其他广义浮点悲伤，我不会感到惊讶。

撇开你的代码不谈，这个 ODE 有大量的术语。它从哪里来的？为了解决这个更具挑战性的问题，您是否可以解决一个更简单的问题？假设您对求解器进行了更改，以便从中得到一些答案。你将如何评估你是否得到了正确的答案？你有确切的解决方案吗？你知道系统是否有任何固定点，如果有，它们的稳定性如何？我问不动点是因为，对于像这样的高次多项式右手边，你很容易有多个根，这会混淆牛顿的方法。

即使在代码中显示的次优版本中，您也在使用牛顿法求解错误的方程。

BDF 方法需要一个隐式方程的解，去掉它可以为一个系统写的分母 $\dot U=F(U)$ 作为

25 \cdot U [I] - 12 \cdot H \cdot F (U [I]) = R = 48 \cdot U [I - 1] - 36 \cdot U [I - 2] + 16 \cdot U [I - 3] - 3 \cdot U [I - 4]

$25·U[I]-12·H·F(U[I]) = R = 48·U[I-1] - 36·U[I-2] + 16·U[I-3] - 3·U[I-4]$ 如果这是一个标量方程，牛顿法解决这个问题

U [I]

$U[I]$ 读作

U_{t m p} = U_{t m p} - \frac{25 \cdot U_{t m p} - 12 \cdot H \cdot F (U_{t m p}) - R}{25 - 12 \cdot H \cdot F^{'} (U_{t m p})}

$U_{\rm tmp} = U_{\rm tmp} - \frac{25·U_{\rm tmp}-12·H·F(U_{\rm tmp}) - R}{25 - 12·H·F'(U_{\rm tmp})}$ 在非标量系统中，导数是雅可比矩阵，计算牛顿步需要求解线性系统。

U_{t m p} = U_{t m p} - L i n S o l v e (25 \cdot I d - 12 \cdot H \cdot F^{'} (U_{t m p}), 25 \cdot U_{t m p} - 12 \cdot H \cdot F (U_{t m p}) - R)

$U_{\rm tmp} = U_{\rm tmp} - {\rm LinSolve}(25·Id - 12·H·F'(U_{\rm tmp}), 25·U_{\rm tmp}-12·H·F(U_{\rm tmp}) - R)$

将系统分离为类似于 Gauß-Seidel 方法的标量方程是可能的，但会破坏牛顿方法的二次收敛性。我认为这不会比仅仅迭代方法方程更好

U [I] = (R + 12.0 \cdot H \cdot F (U [I])) / 25;

$U[I] = (R + 12.0·H·F(U[I])) / 25;$ 直到变化足够小。

作为参考，尽管系数很大，度数很高，但解确实会迅速收敛到一个固定点，

与价值观[2.058626513, 1.0159951728]。这是来自具有相当标准的最小步长控制机制的 Fehlberg45 方法的实现。其他值为

tolerance  steps   final state vector
  1e-04     814    [2.058627514794502,  1.0159952544532254]
  1e-06     917    [2.0586265326480455, 1.0159951744024687]
  1e-08    1258    [2.058626513094132,  1.015995172808779 ]
  1e-10    2336    [2.058626513240043,  1.0159951728206726]

实际上，根据容错要求，有效数字的数量正在增加 2。渐近段上的步长在与h=0.006容错无关的情况下振荡。这是由于系统的刚度和方法是明确的。使用隐式方法，人们会期望在常数段上增加步长。

其它你可能感兴趣的问题

上一篇科学代码使用比较 NaN == NaN 的频率如何？下一篇MATLAB中的匿名函数