godesim

package module
v0.9.3 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Mar 28, 2022 License: BSD-3-Clause Imports: 12 Imported by: 0

README

Go Report Card go.dev reference codecov Awesome

godesim

Simulate complex systems with a simple API.

Wrangle non-linear differential equations while writing maintainable, simple code.

Why Godesim?

ODE solvers seem to fill the niche of simple system solvers in your numerical packages such as scipy's odeint/solve_ivp.

Among these integrators there seems to be room for a solver that offers simulation interactivity such as modifying the differential equations during simulation based on events such as a rocket stage separation.

Installation

Requires Go.

go get github.com/soypat/godesim

Progress

Godesim is in early development and will naturally change as it is used more. The chart below shows some features that are planned or already part of godesim.

Status legend Planned Started Prototype Stable Mature
Legend symbol ✖️ 🏗️ 🐞️ 🚦️ ✅️
Features Status Notes
Non-linear solvers 🚦️ Suite of ODE solvers available.
Non-autonomous support 🚦️ U vector which need not a defined differential equation like X does.
Event driver 🚦️ Eventer interface implemented.
Stiff solver 🚦️ Newton-Raphson algorithm implemented and tested.
Algorithms available and benchmarks
Algorithm Time/Operation Memory/Op Allocations/Op
RK4 1575 ns/op 516 B/op 12 allocs/op
RK5 2351 ns/op 692 B/op 21 allocs/op
RKF45 3229 ns/op 780 B/op 25 allocs/op
Newton-Raphson 8616 ns/op 4292 B/op 92 allocs/op
Dormand-Prince 4365 ns/op 926 B/op 32 allocs/op

Examples

Quadratic Solution
// Declare your rate-of-change functions using state-space symbols
Dtheta := func(s state.State) float64 {
	return s.X("theta-dot")
}

DDtheta := func(s state.State) float64 {
    return 1
}
// Set the Simulation's differential equations and initial values and hit Begin!
sim := godesim.New() // Configurable with Simulation.SetConfig(godesim.Config{...})
sim.SetDiffFromMap(map[state.Symbol]state.Diff {
    "theta":  Dtheta,
    "theta-dot": DDtheta,
})
sim.SetX0FromMap(map[state.Symbol]float64{
    "theta":  0,
    "theta-dot": 0,
})
sim.SetTimespan(0.0, 1.0, 10) // One second simulated
sim.Begin()

The above code solves the following system:

for the domain t=0 to t=1.0 in 10 steps where theta and theta-dot are the X variables. The resulting curve is quadratic as the solution for this equation (for theta and theta-dot equal to zero) is

How to obtain results
// one can then obtain simulation results as float slices 
t := sim.Results("time")
theta := sim.Results("theta")
Other examples

To run an example, navigate to it's directory (under examples) then type go run . in console.

There are three simple examples which have been cooked up and left in _examples directory. I've been having problems running Pixel on my machine so the simulation animations are still under work.

Final notes

Future versions of gonum will have an ODE solver too. Ideally godesim would base it's algorithms on gonum's implementation. See https://github.com/gonum/exp ode package.

Contributing

Pull requests welcome!

This is my first library written for any programming language ever. I'll try to be fast on replying to pull-requests and issues.

Documentation

Overview

Package godesim can be described as a simple interface to solve a first-order system of non-linear differential equations which can be defined as Go code.

Example (Events)

Below is an example of Eventer implementation. `TypicalEventer` implements godesim's Eventer interface. 

type TypicalEventer struct {
	action func(state.State) func(*Simulation) error
	label  string
}
func (ev TypicalEventer) Event(s state.State) func(*Simulation) error { return ev.action(s) }
func (ev TypicalEventer) Label() string { return ev.label }

package main

import (
	"github.com/soypat/godesim"
	"github.com/soypat/godesim/state"
)

type TypicalEventer struct {
	action func(state.State) func(*godesim.Simulation) error
	label  string
}

func (ev TypicalEventer) Event(s state.State) func(*godesim.Simulation) error {
	return ev.action(s)
}

func (ev TypicalEventer) Label() string {
	return ev.label
}

func main() {
	sim := godesim.New()
	sim.SetDiffFromMap(map[state.Symbol]state.Diff{
		"theta": func(s state.State) float64 { return 2 * s.X("theta") },
	})
	sim.SetX0FromMap(map[state.Symbol]float64{
		"theta": 0,
	})

	sim.SetTimespan(0, 10., 10)
	initStepLen := sim.Dt()
	// We halve the step length somewhere along our simulation.
	// this will be the event
	newStepLen := initStepLen * 0.5

	var refiner godesim.Eventer = TypicalEventer{
		label: "refine",
		action: func(s state.State) func(*godesim.Simulation) error {
			if s.Time() >= 3. {
				return godesim.NewStepLength(newStepLen)
			}
			return nil
		},
	}
	sim.AddEventHandlers(refiner)
	sim.Solver = godesim.NewtonRaphsonSolver
	sim.Begin()
}

Output:
Example (Implicit)

Solve a stiff equation problem using the Newton-Raphson method. Equation being solved taken from https://en.wikipedia.org/wiki/Stiff_equation Do note that the accuracy for stiff problems is reduced greatly for conventional methods. 

y'(t) = -15*y(t)
solution: y(t) = exp(-15*t)

package main

import (
	"fmt"
	"math"

	"github.com/soypat/godesim"
	"github.com/soypat/godesim/state"
)

func main() {
	sim := godesim.New()
	tau := -15.
	solution := func(x []float64) []float64 {
		sol := make([]float64, len(x))
		for i := range x {
			sol[i] = math.Exp(tau * x[i])
		}
		return sol
	}
	sim.SetDiffFromMap(map[state.Symbol]state.Diff{
		"y": func(s state.State) float64 {
			return tau * s.X("y")
		},
	})
	sim.SetX0FromMap(map[state.Symbol]float64{
		"y": 1,
	})
	sim.SetTimespan(0.0, 0.5, 15)
	sim.Begin()
	fmt.Printf("domain  :%0.3f:\nresult  :%0.3f\nsolution:%0.3f", sim.Results("time"), sim.Results("y"), solution(sim.Results("time")))
}

Output:

domain  :[0.000 0.033 0.067 0.100 0.133 0.167 0.200 0.233 0.267 0.300 0.333 0.367 0.400 0.433 0.467 0.500]:
result  :[1.000 0.607 0.368 0.223 0.136 0.082 0.050 0.030 0.018 0.011 0.007 0.004 0.002 0.002 0.001 0.001]
solution:[1.000 0.607 0.368 0.223 0.135 0.082 0.050 0.030 0.018 0.011 0.007 0.004 0.002 0.002 0.001 0.001]
Example (Quadratic)

Solves a simple system of equations of the form 

Dtheta     = theta_dot
Dtheta_dot = 1

package main

import (
	"fmt"

	"github.com/soypat/godesim"
	"github.com/soypat/godesim/state"
)

func main() {
	sim := godesim.New()
	sim.SetDiffFromMap(map[state.Symbol]state.Diff{
		"theta": func(s state.State) float64 {
			return s.X("theta-dot")
		},
		"theta-dot": func(s state.State) float64 {
			return 1
		},
	})
	sim.SetX0FromMap(map[state.Symbol]float64{
		"theta":     0,
		"theta-dot": 0,
	})
	sim.SetTimespan(0.0, 1.0, 10)
	sim.Begin()
	fmt.Printf("%0.3f:\n%0.3f", sim.Results("time"), sim.Results("theta"))
}

Output:

[0.000 0.100 0.200 0.300 0.400 0.500 0.600 0.700 0.800 0.900 1.000]:
[0.000 0.005 0.020 0.045 0.080 0.125 0.180 0.245 0.320 0.405 0.500]

Index

Examples

Constants

This section is empty.

Variables

View Source 
var (
	// ErrorRemove should be returned by a simulation modifier
	// if the Eventer is to be removed from Event handlers
	ErrorRemove error = fmt.Errorf("remove this Eventer")
)

Functions

func DiffChangeFromMap 

func DiffChangeFromMap(newDiff map[state.Symbol]func(state.State) float64) func(*Simulation) error

DiffChangeFromMap Event handler. Takes new state variable (X) functions and applies them

func DirectIntegrationSolver added in v0.9.1 

func DirectIntegrationSolver(sim *Simulation) []state.State

DirectIntegrationSolver performs naive integration of ODEs. Should only be used to compare with other methods. Has the advantage of only performing one differentiation per step. 

y_{n+1} = y_{n} + (dy/dt_{n} + dy/dt_{n-1})*step/2

func DormandPrinceSolver added in v0.8.1 

func DormandPrinceSolver(sim *Simulation) []state.State

DormandPrinceSolver. Very similar to RKF45. Used by Matlab in ode45 solver and Simulink's system solver by default.

To enable adaptive stepping, Config.Algorithm.Step Min/Max values must be set and a Config.Error.Min must be specified in configuration.

func EndSimulation 

func EndSimulation(sim *Simulation) error

EndSimulation Event handler. Ends simulation

func EventDone 

func EventDone(sim *Simulation) error

EventDone is the uneventful event. Changes absolutely nothing and is removed from Event handler list.

func NewStepLength 

func NewStepLength(h float64) func(*Simulation) error

NewStepLength Event handler. Sets the new minimum step length

func NewtonRaphsonSolver added in v0.7.0 

func NewtonRaphsonSolver(sim *Simulation) []state.State

NewtonRaphsonSolver is an implicit solver which may calculate the jacobian several times on each algorithm step.

sim.Algorithm.Error.Max should be set to a value above 0 for good run

func RK4Solver 

func RK4Solver(sim *Simulation) []state.State

RK4Solver Integrates simulation state for next timesteps using 4th order Runge-Kutta multivariable algorithm

func RKF45Solver 

func RKF45Solver(sim *Simulation) []state.State

RKF45Solver Runge-Kutta-Fehlberg of Orders 4 and 5 solver

To enable adaptive stepping, Config.Algorithm.Step Min/Max values must be set and a Config.Error.Min must be specified in configuration.

func RKF78Solver added in v0.9.2 

func RKF78Solver(sim *Simulation) []state.State

func StateDiff 

func StateDiff(F state.Diffs, S state.State) state.State

StateDiff obtain Diffs results without modifying State Returns state evolution (result of applying Diffs functions to S)

Types

type Config 

type Config struct {
	// Domain is symbol name for step variable. Default is `time`
	Domain    state.Symbol `yaml:"domain"`
	Log       LoggerOptions
	Behaviour struct {
		StepDelay time.Duration `yaml:"delay"`
	} `yaml:"behaviour"`
	Algorithm struct {
		// Number of algorithm steps. Different to simulation Timespan.Len()
		Steps int `yaml:"steps"`
		// Step limits for adaptive algorithms
		Step struct {
			Max float64 `yaml:"max"`
			Min float64 `yaml:"min"`
		} `yaml:"step"`
		Error struct {
			// Sets max error before proceeding with adaptive iteration
			// Step.Min should override this
			Max float64 `yaml:"max"`
		} `yaml:"error"`

		// Newton-Raphson method's convergence can benefit from a relaxationFactor between 0 and 0.5
		RelaxationFactor float64 `yaml:"relaxation"`
		// Newton-Raphson Method requires multiple sub-iterations to converge. On each
		// iteration the Jacobian is calculated, which is an expensive operation.
		// A good number may be between 10 and 100.
		IterationMax int `yaml:"iterations"`
	} `yaml:"algorithm"`
	Symbols struct {
		// Sorts symbols for consistent logging and testing
		NoOrdering bool `yaml:"no_ordering"`
	} `yaml:"symbols"`
}

Config modifies Simulation behaviour/output. Set with simulation.SetConfig method

func DefaultConfig added in v0.7.3 

func DefaultConfig() Config

DefaultConfig returns configuration set for all new simulations by New()

Domain is the integration variable. "time" is default value 

simulation.Domain

Solver used is fourth order Runge-Kutta multivariable integration. 

simulation.Solver

How many solver steps are run between Timespan steps. Set to 1 

simulation.Algorithm.Steps

type Eventer 

type Eventer interface {
	// Behaviour aspect of event. changing simulation
	//
	// A nil func(*Simulation) error corresponds to no event taking place (just skips it)
	Event(state.State) func(*Simulation) error
	// For identification
	Label() string
}

Eventer specifies an event to be applied to simulation

Events which return error when applied to simulation will create a special event with the error message in the label and log the information.

type Logger 

type Logger struct {
	Output io.Writer
	// contains filtered or unexported fields
}

Logger accumulates messages during simulation run and writes them to Output once simulation finishes.

func (*Logger) Logf 

func (log *Logger) Logf(format string, a ...interface{})

Logf formats message to simulation logger. Messages are printed when simulation finishes. This is a rudimentary implementation of a logger.

type LoggerOptions added in v0.7.0 

type LoggerOptions struct {
	Results struct {
		AllStates bool
		FormatLen int    `yaml:"format_len"`
		Separator string `yaml:"separator"`
		// Defines X in formatter %a.Xg for floating point values.
		// A value of -1 specifies default precision
		Precision int `yaml:"prec"`
	} `yaml:"results"`
}

LoggerOptions for now permits user to output formatted values to an io.Writer.

Example of usage for csv generation: 

logcfg := godesim.LoggerOptions{}
logcfg.Results.Separator = ","
logcfg.Results.AllStates = true
logcfg.Results.FormatLen = 6
cfg := godesim.DefaultConfig()
cfg.Log = logcfg

Finally, set simulation config with sim.SetConfig(cfg) method and set io.Writer in sim.Logger.Output

type Simulation 

type Simulation struct {
	Timespan
	Logger Logger
	State  state.State

	Solver func(sim *Simulation) []state.State

	Diffs state.Diffs

	Config
	// contains filtered or unexported fields
}

Simulation contains dynamics of system and stores simulation results.

Defines an object that can solve a non-autonomous, non-linear system of differential equations

func New 

func New() *Simulation

New creates blank simulation. To run a simulation one must also set the domain with NewTimespan() and create the differential equation system with SetX0FromMap() and SetDiffFromMap(). Examples can be found at https://pkg.go.dev/github.com/soypat/godesim.

The Default solver is RK4Solver. Other default values are set by DefaultConfig()

func (*Simulation) AddEventHandlers 

func (sim *Simulation) AddEventHandlers(evhand ...Eventer)

AddEventHandlers add event handlers to simulation.

Events which return errors will create a special event with an error message for Label(). If one wishes to stop simulation execution one can call panic() in an event.

func (*Simulation) Begin 

func (sim *Simulation) Begin()

Begin starts simulation

Unrecoverable errors will panic. Warnings may be printed.

func (*Simulation) CurrentTime 

func (sim *Simulation) CurrentTime() float64

CurrentTime obtain simulation step variable

func (*Simulation) Events 

func (sim *Simulation) Events() []struct {
	Label string
	State state.State
}

Events Returns a copy of all simulation events

func (*Simulation) ForEachState added in v0.9.0 

func (sim *Simulation) ForEachState(f func(i int, s state.State))

ForEachState calls f on all result states. Simulation must have been run beforehand.

func (*Simulation) IsRunning 

func (sim *Simulation) IsRunning() bool

IsRunning returns true if Simulation has yet not run last iteration

func (*Simulation) Results 

func (sim *Simulation) Results(sym state.Symbol) []float64

Results get vector of simulation results for given symbol (X or U)

Special case is the Simulation.Domain (default "time") symbol.

func (*Simulation) SetConfig 

func (sim *Simulation) SetConfig(cfg Config) *Simulation

SetConfig Set configuration to modify default Simulation values

func (*Simulation) SetDiffFromMap 

func (sim *Simulation) SetDiffFromMap(m map[state.Symbol]state.Diff)

SetDiffFromMap Sets the ODE equations (change in X) with a pre-built map

i.e. theta(t) = 0.5 * t^2 

sim.SetDiffFromMap(map[state.Symbol]state.Change{
	"theta":  func(s state.State) float64 {
		return s.X("Dtheta")
	},
	"Dtheta": func(s state.State) float64 {
		return 1
	},
})

func (*Simulation) SetInputFromMap 

func (sim *Simulation) SetInputFromMap(m map[state.Symbol]state.Input)

SetInputFromMap Sets Input (U) functions with pre-built map

func (*Simulation) SetX0FromMap 

func (sim *Simulation) SetX0FromMap(m map[state.Symbol]float64)

SetX0FromMap sets simulation's initial X values from a Symbol map

func (*Simulation) States added in v0.9.0 

func (sim *Simulation) States() (states []state.State)

StatesCopy returns a copy of all result states. Simulation must have been run beforehand.

type Timespan 

type Timespan struct {
	// contains filtered or unexported fields
}

Timespan represents an iterable vector of evenly spaced time points. Does not store state information on steps done.

func (Timespan) Dt 

func (ts Timespan) Dt() float64

Dt Obtains the step length of simulation

func (Timespan) End 

func (ts Timespan) End() float64

End returns greater limit of Timespan

func (Timespan) Len 

func (ts Timespan) Len() int

Len how many iterations expected for RK4

func (*Timespan) SetTimespan 

func (ts *Timespan) SetTimespan(Start, End float64, Steps int)

SetTimespan Set time domain (step domain) for simulation. Step size is given by: 

dt = (End - Start) / float64(Steps)

since Steps is the amount of points to "solve".

Source Files

View all Source files

Directories

Path Synopsis
_examples
doublePendulum
equations from https://www.math24.net/double-pendulum/
 equations from https://www.math24.net/double-pendulum/
n-body
simplePendulum

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.