Quick start to ODE

The following code shows how to obtain the numerical solution of a simple one-dimensional ODE. Note the configuration of the stepper type showing the flexibility of the library. We solve the simple ODE x' = 3/(2t^2) + x/(2t) with initial condition x(1) = 0 and step size 0.1.

dt := 0.1. 
system := PMExplicitSystem block: [:x :t | 3.0 / (2.0 * t * t) + (x / (2.0 * t))]. 
stepper := PMRungeKuttaStepper onSystem: system. 
solver := (PMExplicitSolver new) stepper: stepper; system: system; dt: dt. 
solver solve: system startState: 0 startTime: 1 endTime: 10

Every time, you need to solve an ODE, you will need to setup three different object: a system that represent an or a set of ODEs, a stepper that deals with the method that will be used to solve your ODE and finally a solver.

With the PMExplicitSystem block you can represent as a function of 2 arguments (t and x), the rhs (right hand side) of the differential equation: dx/dt = f(x,t).

The analytic solution is x(t) = sqrt(t) - 1/t

Let's try to do plot the behavior of an compartimental model in epidemiology like SIR:

dt := 1.0.
beta := 0.01.
gamma := 0.1.
system := PMExplicitSystem block: [ :x :t| {
	(beta negated) * (x at: 1) * (x at: 2).
	(beta * (x at: 1) * (x at: 2)) - (gamma * (x at: 2)).
	gamma * (x at: 2).
	}].
stepper := PMRungeKuttaStepper onSystem: system.
solver := (PMExplicitSolver new) stepper: stepper; system: system; dt: dt.
state := #(99 1 0).
values := (0.0 to: 75.0 by: dt) collect: [ :t| state := stepper doStep: state 
                                                          time: t stepSize: dt. t@(state at:2)].

plt := RSChart new.
x := values collect: #x.
y := values collect: #y.

p := RSLinePlot new x: x y: y.
plt addPlot: p.

plt addDecoration: RSHorizontalTick new.
plt addDecoration: RSVerticalTick new.
plt open

We use Roassal3, to plot the S curve.

For the three curves: S (blue), I (red) and R (green)

|solver state system dt beta gamma values stepper grapher dataSet|
dt := 1.0.
beta := 0.01.
gamma := 0.1.
system := PMExplicitSystem block: [ :x :t| {
	(beta negated) * (x at: 1) * (x at: 2).
	(beta * (x at: 1) * (x at: 2)) - (gamma * (x at: 2)).
	gamma * (x at: 2).
	}].
stepper := PMRungeKuttaStepper onSystem: system.
solver := (PMAB4Solver new) stepper: stepper; system: system; dt: dt.
state := #(99 1 0).
values := (0.0 to: 60.0 by: dt) collect: [ :t| state := stepper doStep: state 
                                                          time: t stepSize: dt. {t. state at:1. state at:2. state at:3}].
grapher := RTGrapher new.

dataSet := RTData new.
dataSet noDot.
dataSet points: values;
	x: [:v| v at:1];
	y: [:v| v at:2].
dataSet connectColor: (Color green).
grapher add: dataSet.

dataSet := RTData new.
dataSet noDot.
dataSet points: values;
	x: [:v| v at:1];
	y: [:v| v at:3].
dataSet connectColor: (Color red).
grapher add: dataSet.

dataSet := RTData new.
dataSet noDot.
dataSet points: values;
	x: [:v| v at:1];
	y: [:v| v at:4].
dataSet connectColor: (Color blue).
grapher add: dataSet.

grapher

Let's try to implement the Lorenz attractor as the a set of 3 differential equations, that we will solved with a Runge-Kutta solver:

sigma := 10.0.
r := 28.
b := 8.0/3.0.
dt := 0.01.
system := PMExplicitSystem
block: [:x :t | | c |
   c:= Array new: 3.
   c at: 1 put: sigma * ((x at: 2) - (x at: 1)).
   c at: 2 put: r * (x at: 1) - (x at: 2) - ((x at: 1) * (x at: 3)).
   c at: 3 put: (b negated * (x at: 3) + ((x at: 1) * (x at: 2))).
   c].
stepper := PMRungeKuttaStepper onSystem: system.
solver := (PMExplicitSolver new) stepper: stepper; system: system; dt: dt.
state := #(10.0 10.0 10.0).

values := (0.0 to: 100.0 by: dt) collect:[:t |
   state := stepper doStep: state time: t stepSize: dt.
   state].

grapher := RTGrapher new.
dataSet := RTData new.
dataSet dotShape ellipse size:1; color: (Color red alpha: 0.5).
dataSet
	points: values;
	x: [:v | v at: 2];
	y: [:v | v at: 3].
grapher add: dataSet.
grapher.

As there are different case what value to write for system, stepper and solver here you have a hint:

(Pay attentions: for first 3 methods you should write your order of method like thisAB3Stepper. )

Methods	System	Stepper	Solver	Comment
Adams - Bashforth method of order n=2,3,4	ExplicitSystem	ABnStepper	ABnSolver
Adams - Moulton method of order n=3,4	ImplicitSystem	AMnStepper	AMnSolver
Backward differentiation formulas method of order n=2,3,4	ImplicitSystem	BDFnStepper	BDFnSolver
Beckward Euler method	ImplicitSystem	ImplicitStepper	ImplicitSolver
Euler method	ExplicitSystem	ExplicitStepper	ExplicitSolver
Heun method	ExplicitSystem	HeunStepper	ExplicitSolver	Heun's method may refer to the improved or modified Euler's method (that is, the explicit trapezoidal rule), or a similar two-stage Runge-Kutta method.
Implicit midpoint method	ImplicitSystem	ImplicitMidpointStepper	ImplicitMidpointSolver
Midpoint method	ExplicitSystem	MidpointStepper	ExplicitSolver
Runge-Kutta method (order 4)	ExplicitSystem	RungeKuttaStepper	ExplicitSolver
Trapezoidal rule	ImplicitSystem	TrapezoidStepper	ImplicitSolver	It is Adams - Moulton method of order 2