Demonstration of class ODESolver usage: A linear 2nd order differential system of equations

We present a program to solve by a numerical scheme a system of first-order linear differential equations of the form 

                  A2y'' + A1y' + A0y = f
where the matrices A0 and A1 and A2 are 
        | 2  -1 |         | -2  0 |       | 2  0 |
   A0 = |       |,   A1 = |       |, A2 = |      |
        | -1  2 |         |  0  1 |       | 0  1 |
and 
       |  (2-2π2)sin(πt) - 2πcos(πt) - cos(2πt) |
   f = |                                        |.
       | (2-4π2)cos(2πt) - 2πsin(2πt) - sin(πt) |
With adapted initial conditions, the solution is given by 
        y1(t) = sin(πt), y2 = cos(2πt)
  • After initializing the global variables theFinalTime and theTimeStep, we declare an instance of class ODESolver, where we have chosen the Newmark scheme with default parameters, and set the number of equations to 2. 

       ODESolver ode(NEWMARK,theTimeStep,theFinalTime,2);

  • Once the instance created, we declare the matrices (as dense matrices) and send them to the class. Note that since their pointers are stored, their contents can be modified thrghout the time stepping. Next we declare the necessary vectors: the solution and the right-hand side vector. We give the initial condition using the member function setInitial, the right-hand side using setRHS and, since we deal with a two-step method, an initial right-hand side. Note that this is not mandatory but improves the accuracy. 

       DMatrix<double> A0(2,2), A1(2,2), A2(2,2);
   ode.setMatrices(A0,A1,A2);
   Vect<double> u(2), v(2);
   u(1) = 0; u(2) = 1;
   v(1) = OFELI_PI; v(2) = 0;
   ode.setInitial(u,v);
   ode.setRHS("(2-2*pi^2)*sin(pi*t)-2*pi*cos(pi*t)-cos(2*pi*t)");
   ode.setRHS("(2-4*pi^2)*cos(2*pi*t)-2*pi*sin(2*pi*t)-sin(pi*t)");

  • Since the system of equations as well as data are given, in this example, by numerical values (rather than by algebraic expressions), we explicitly use a loop over time steps, where at each time step, the matrices, as well as the right-hand side are updated. Note that this is necessary even if the matrices are constant, since they are fatorized each time step. 

       TimeLoop {
      A0(1,1) =  2.; A0(1,2) = -1.; A0(2,1) = -1.; A0(2,2) =  2.;
      A1(1,1) = -2.; A1(1,2) =  0.; A1(2,1) =  0.; A1(2,2) =  1.;
      A2(1,1) =  2.; A2(1,2) =  0.; A2(2,1) =  0.; A2(2,2) =  1.;
      ode.runOneTimeStep();
   }

  • The solution is stored each time step in the vector u. We can in addition compute the error: 

       cout << ode << endl;
   cout << "Error: " << fabs(u(1)-sin(OFELI_PI*theFinalTime)) << ", "
        << fabs(u(2)-cos(2*OFELI_PI*theFinalTime)) << endl;
   return 0;
 }