Skip to content

How does ODE package work?

Serge Stinckwich edited this page Jun 26, 2018 · 1 revision

To understand how you can work with the ODE package, let's look at the example with Adams-Bashforth second order method. Read about Adams-Bashforth method and see what the algorithm we will use. Let's start. First of all you need to understand how this method works. If you want to solve some ODE using this method you need to do next steps:

  1. Check whether t1 is equal to lastT, if not then
  2. Find the first state x1 using first order method from Adams-Bashforth family of methods and initial data: initial state x0, initial time t0 and step size dt. Check whether t2 is equal to lastT, if not then
  3. Find the second state x2 using second order Adams-Bashforth method and previous datas: x1, x0, t1, t0, dt. Check whether t3 is equal to lastT, if not then
  4. Find the n next states x3,...,xn using second order Adams-Bashforth method and previous datas. Here n is equal to (t-t0)/dt. If on some (k+1) step you have tk<lastT<t(k+1) then find x(k+1)_ using k+1 order Adams-Bashforth method but with new_dt=lastT-tk. If on some k step tk=lastT then find xk and it's all.

Now, when you know how does the method works, it's a time to look on the AB2Solver. The solve method do all previous steps.

solve: aSystem startState: initialState startTime: initialTime endTime: endTime
|prevState statesPair|
self system: aSystem.
self stepper: ((self firstStepperClass) onSystem: self system).
state := initialState.
lastTime:=initialTime.

"announce initial conditions"
self announceState: state time: initialTime.

(lastTime+dt > endTime and: lastTime < endTime)
ifTrue:[
    state :=self lastStepState: state endTime: endTime].

lastTime+dt<= endTime 
ifTrue:[
    prevState:= initialState.

    state := self firstStepStartTime: initialTime.
    "announce first step"
    self announceState: state time: (initialTime + dt).

    self stepper: ((self stepperClass) onSystem: self system).

    "step until the end"
    statesPair := self mainStepsPrevState: prevState startTime: ( initialTime + dt) endTime: endTime.
    state:=statesPair second.
    prevState:= statesPair first.

    "sanity check"
    self assert: [(lastTime between: initialTime and: endTime) 
        or: [lastTime between: endTime and: initialTime]].

    "take another step if needed"
    state := self lastStepPrevState: prevState endTime: endTime.].

^ state

In previous code-block we used next methods:

firstStepStartTime: t
  state := stepper doStep: state time: t stepSize: self dt. 
  self announceState: state time: t + self dt.
  lastTime := t + self dt.
  ^ state

which sent a message to MidpointStepper. This method does the second step.

mainStepsPrevState: prevState startTime: initialTime endTime: endTime
|previousState|
previousState := prevState. 
"don't go to end time to avoid overrunning"
(initialTime to: endTime - self dt  by: self dt) do:
[:time | | tempState|
    tempState := state.
    state := stepper
        doStep: state 
        prevState: previousState
        time:  time 
        stepSize: self dt.
    previousState := tempState.
    "announce step results"
    self announceState: state time: time + self dt.
    lastTime := time + self dt].        
^ {previousState . state}

This method does the third and forth steps. It sents the messages to the method doStep: prevState: from AB2Stepper

doStep: aState prevState: prevState time: t
    self stepSize isNil
     ifTrue: [ self error: 'step size required by stepper' ].
    ^ (self stepSize / 2) * (3 * (system x: aState t: t) - (system x: prevState t: t - self stepSize)) + aState


lastStepPrevState: prevState endTime: endTime
"catch partial or full step at end"
    (lastTime equalsTo: endTime ) ifFalse:
    [state := stepper 
    lastStep: state 
    prevState: prevState
    time: lastTime 
    stepSize: endTime - lastTime
    deltaT: dt.
    self announceState: state time: endTime].           
^ state

This method works in this case tk<lastT<t(k+1). It also sents a message to the method lastStep: prevState: time: deltaT: from AB2Stepper

lastStep: aState prevState: prevState time: t deltaT: incrementOfTime
    self stepSize isNil
    ifTrue: [ self error: 'step size required by stepper' ].
    ^ self stepSize / 2 * (3 * (system x: aState t: t) - (system x: prevState t: t - incrementOfTime)) + aState