Design of an ODE Solver Environment

1. Design of an ODE Solver Environment

2. System of ODEs Consider the system of ODEs: Typically solved by forward Euler or 4th order Runge-Kutta

3. Traditional Approach Traditionally solved by library routines: SUBROUTINE RK4 (Y,T,F,WORK1,N,TSTEP,TOL1,TOL2,…) Y: current value of y T: current value for time t F: external function defining f WORK1: work array N: length of array TSTEP: time step value TOL1,TOL2:algorithmic parameters

4. Model Problem Define and you get:

5. Traditional Interface That is, the function F has the following interface: SUBROUTINE F (YDOT,Y,T,C1,C2,C3,C4,OMEGA) NB! Problem-dependent interface!! Solvers like RK4 wants: SUBROUTINE F (YDOT,Y,T) (Bad) Solution: COMMON block (global data)

6. The Object-Oriented Way • User program can access base class (solver unknown) • Class hierarchy of solvers • Embedding of problem-specific parameters • Uniform interface for function f

7. ForwardEuler RungeKutta2 ODEProblem ODESolver RungeKutta4 RungeKutta4A Oscillator ... ... Suggested Design

8. Class ODESolver class ODESolver { protected: // members only visible in subclasses ODEProblem* eqdef; // definition of the ODE in user's class public: // members visible also outside the class ODESolver (ODEProblem* eqdef_) { eqdef = eqdef_; } virtual ~ODESolver () {} // always needed, does nothing here... virtual void init(); // initialize solver data structures virtual void advance (Vec(real)& y, real& t, real& dt); };

9. Class ODEProblem class ODEProblem { protected: ODESolver* solver; // some ODE solver Vec(real) y, y0; // solution (y) and initial condition (y0) real t, dt, T; // time loop parameters public: ODEProblem () {} virtual ~ODEProblem (); virtual void timeLoop (); virtual void equation (Vec(real)& f, const Vec(real)& y, real t); virtual int size (); // no of equations in the ODE system virtual void scan (); virtual void print (Os os); };

10. The Time Loop void ODEProblem:: timeLoop () { Os outfile (aform("%s.y",casename.c_str()), NEWFILE); t = 0; y = y0; outfile << t << " "; y.print(outfile); outfile << '\n'; while (t <= T) { solver->advance (y, t, dt); // update y, t, and dt outfile << t << " "; y.print(outfile); outfile << '\n'; } }

11. Class Oscillator class Oscillator : public ODEProblem { protected: real c1,c2,c3,c4,omega; // problem dependent paramters public: Oscillator () {} virtual void equation (Vec(real)& f, const Vec(real)& y, real t); virtual int size () { return 2; } // 2x2 system of ODEs virtual void scan (); virtual void print (Os os); };

12. Defining the ODE void Oscillator::equation (Vec(real)& f, const Vec(real)& y_, real t_) { f(1) = y_(2); f(2) = -c1*(y_(2)+c2*y_(2)*abs(y_(2))) - c3*(y_(1)+c4*pow3(y_(1))) + sin(omega*t_); }

13. Class ForwardEuler class ForwardEuler : public ODESolver { Vec(real) scratch1; // needed in the algorithm public: ForwardEuler (ODEProblem* eqdef_); virtual void init (); // for allocating scratch1 virtual void advance (Vec(real)& y, real& t, real& dt); };

14. Stepping the Solver void ForwardEuler:: advance (Vec(real)& y, real& t, real& dt) { eqdef->equation (scratch1, y, t); // evaluate scratch1 (as f) const int n = y.size(); for (int i = 1; i <= n; i++) y(i) += dt * scratch1(i); t += dt; }

15. Parameter Classes class ODESolver_prm { public: String method; // name of subclass in ODESolver hierarchy ODEProblem* problem; // pointer to user's problem class ODESolver* create (); // create correct subclass of ODESolver }; ODESolver* ODESolver_prm:: create () { ODESolver* ptr = NULL; if (method == "ForwardEuler") ptr = new ForwardEuler (problem); else if (method == "RungeKutta4") ptr = new RungeKutta4 (problem); else errorFP("ODESolver_prm::create", "Method \"%s\" is not available",method.c_str()); return ptr; }

16. Problem Input void ODEProblem:: scan () { const int n = size(); // call size in actual subclass y.redim(n); y0.redim(n); s_o << "Give " << n << " initial conditions: "; y0.scan(s_i); s_o << "Give time step: "; s_i >> dt; s_o << "Give final time T: "; s_i >> T; ODESolver_prm solver_prm; s_o << "Give name of ODE solver: "; s_i >> solver_prm.method; solver_prm.problem = this; solver = solver_prm.create(); solver->init(); // more reading in user's subclass }

17. Problem Input void Oscillator:: scan () { // first we need to do everything that ODEProblem::scan does: ODEProblem::scan(); // additional reading here: s_o << "Give c1, c2, c3, c4, and omega: "; s_i >> c1 >> c2 >> c3 >> c4 >> omega; print(s_o); // convenient check for the user }

18. The Main Program int main (int argc, const char* argv[]) { initDiffpack (argc, argv); Oscillator problem; problem.scan(); // read input data and initialize problem.timeLoop(); // solve problem return 0; }

19. What Has Been Gained? • Some “disturbance” of basic numerics • Flexible library • Faster to define problem (f) than hardcoding solver • Simple example, ideas transfer to very complicated scenarios

20. Exercise: Extending the ODE Environment Use the ODE environment to solve the problem u’(t) = -u(t) u(0) = 1 The analytical solution to this problem is u(t)=e-t.