- 54 Views
- Uploaded on
- Presentation posted in: General

Design of an ODE Solver Environment

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Design of an ODE Solver Environment

System of ODEs

Consider the system of ODEs:

Typically solved by forward Euler or 4th order Runge-Kutta

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

Model Problem

Define

and you get:

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)

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

ForwardEuler

RungeKutta2

ODEProblem

ODESolver

RungeKutta4

RungeKutta4A

Oscillator

...

...

Suggested Design

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);

};

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);

};

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';

}

}

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);

};

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_);

}

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);

};

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;

}

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;

}

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

}

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

}

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;

}

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

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.