3
$\begingroup$

I'm trying to solve a system of ODEs using the BDF order 4 method. I find the first 3 points using RK4, then for the implicit part of the BDF, I use Newton-Raphson iteration. Unfortunately my solution blows up but I can't figure what I'm doing wrong. I assume it's something to do with the Newton iteration since I'm not too confident about that part. The ODE is of the form

$y_1' = 40·y_1 -2·y_2^2 +40·y_2^2 -100·y_1·y_2^2+...$ about 30 more terms

$y_2' = 2·y_1 +40·y_2 +...$ about 30 more terms

ODEs

Here's my code so far

#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;
}
Improve this question
$\endgroup$
6
  • $\begingroup$ Could you please edit your post and write the ODE you're trying to integrate? @MichealB $\endgroup$
    – VoB
    Commented Dec 29, 2020 at 17:03
  • $\begingroup$ Thanks for the reply. I've edited the post. The ODE is a system of 2 equations in y1 and y2 with about 30 terms each, hence I rather not type out the whole thing $\endgroup$
    – Michael B
    Commented Dec 29, 2020 at 17:22
  • $\begingroup$ What do you mean about more thirty terms? I think still you did not addressed @VoB's comment yet. Please use LaTeX and write the complete set of equations no matter how many terms you have. $\endgroup$
    – Mithridates the Great
    Commented Dec 29, 2020 at 18:44
  • $\begingroup$ Well they're just polynomial terms in y1 and y2. I've added an image of the whole equations to the OP. I was under the impression that new accounts couldn't post pictures. Sorry. $\endgroup$
    – Michael B
    Commented Dec 29, 2020 at 18:52
  • 1
    $\begingroup$ Please clean your code up so that it is either C++ or C. A strong hint that this is not the case at the moment is the use of C headers like math.h instead of the C++ header cmath. $\endgroup$
    – Lutz Lehmann
    Commented Jan 2, 2021 at 11:14

2 Answers 2

Reset to default
3
$\begingroup$

A few things jump out at me from your code as potential problems.

  1. You seem to be using Newton's method separately for the $y_1$ and $y_2$ variables. This is not the same as using Newton's method for a nonlinear system involving both $y_1$ and $y_2$. For the full Newton method, you'll be solving a 2x2 linear system at every step. You'll have to calculate not just $\partial F/\partial y_1$ and $\partial S/\partial y_2$ but also $\partial F/\partial y_2$ and $\partial S/\partial y_1$. What you've written down is a coordinate-descent type approach, which can work for some problems, but which nonetheless is not Newton's method.
  2. Newton's method is not guaranteed to converge for any initial data that you throw at it. If the initial guess is very far outside the quadratic convergence basin, the iterates could fly off to infinity. In general, you need a globalization strategy like a line search or trust region method.
  3. The coefficients in your problem are large enough that, for the time interval you're solving over, I wouldn't be surprised if you encounter overflow or other generalized floating-point sadness.

Leaving your code aside, this ODE has a huge number of terms. Where did it come from? Is there a simpler problem you could solve instead in order to work your way up to this more challenging problem? Let's say that you make changes to your solver so that you do get some answer out of it. How will you assess whether you got the right answer? Do you have an exact solution? Do you know if the system has any fixed points and if so what their stability is like? I ask about fixed points because, for a high-degree polynomial right-hand side like this, you could easily have multiple roots that will confound Newton's method.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for the reply. You're definitely right about the Newton method implementation. I found a version that works for non-linear systems using the Jacobian, which is probably what you've suggested. I've used a matlab solver and the solution seems to be decently well behaved (the curves appear sinusoidal in nature), so I hoped convergence wouldn't be a big issue. It's possible that I might have to look into what you suggested. $\endgroup$
    – Michael B
    Commented Dec 29, 2020 at 21:21
  • $\begingroup$ Regarding the ODE, it's from a 'challenge' problem set I have in one of my math classes. I know the solution value at the end of my time step, and I'm supposed to match it with a given accuracy. The solution does have fixed points, but I've read that Newton's method is much more appropriate when solving stiff systems, hence why I tried implementing it. $\endgroup$
    – Michael B
    Commented Dec 29, 2020 at 21:23
  • $\begingroup$ Given that your system is only made of 2 variables, you can simply use an explicit time integrator with adaptive time stepping (these are relatively easy to implement), even if the system is stiff. That might save you a lot of trouble with the implicit terms. The explicit time integration may however be more costly if you are only searching for a long term / steady-state behaviour. You can try solving your problem with existing integrator libraries (Boost or maybe some code found on Github) to see what the solution should look like, and test implicit or explicit integrators. $\endgroup$
    – Laurent90
    Commented Dec 30, 2020 at 13:02
  • $\begingroup$ Unfortunately I'm not allowed to use third party libraries for this problem. After using the Jacobian method I'm still unable to arrive at the correct answer. Since the solution is sinusoidal in nature, do you think it would be possible to just use a high order explicit method instead? $\endgroup$
    – Michael B
    Commented Dec 30, 2020 at 14:01
  • $\begingroup$ In what way is the Newton method failing to give the correct answer? Is it exploding, or is it producing a trajectory that doesn't end where it should? It's worth trying even a low-order explicit method like Euler. Since you know where the system should end up, you can compute the eigenvalues of the linearization of the RHS around the starting and final states to get an idea of what sort of timestep you need. $\endgroup$
    – Daniel Shapero
    Commented Dec 30, 2020 at 18:07
1
$\begingroup$

You are solving the wrong equations with the Newton method, even in the sub-optimal version presented in the code.

The BDF method requires the solution of an implicit equation, removing denominators it can be written for a system $\dot U=F(U)$ as $$ 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] $$ If this were for a scalar equation, the Newton method to solve this for $U[I]$ would read as $$ 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})} $$ In a non-scalar system the derivative is a Jacobi matrix and computing the Newton step requires solving a linear system. $$ 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) $$

Separation of the system into scalar equations analogous to the Gauß-Seidel method is possible, but will destroy the quadratic convergence of the Newton method. I think this would be no better than just iterating the method equation $$ U[I] = (R + 12.0·H·F(U[I])) / 25; $$ until the changes are sufficiently small.

For reference, despite large coefficients and high degrees, the solution indeed rapidly converges to a stationary point,

enter image description here

with values [2.058626513, 1.0159951728]. This is from an implementation of the Fehlberg45 method with a rather standard minimal step size control mechanism. Other values are

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]

where indeed the number of valid digits is progressing by 2 as the error tolerance demanded. The step size on the asymptotic segment oscillated around h=0.006 independent of the error tolerance. This is due to the stiffness of the system and that the method is explicit. With an implicit method one would expect increasing step sizes on the constant segment.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.