Fb1_ODEdef container for ordinary differential equation definitions
Part of: miSCellaneous
To be used to define ODE systems that can then be audified with Fb1_ODE. See Fb1_ODE help for general information about Fb1 ODE integrator UGens.
HISTORY AND CREDITS: Big credit to David Pirrò from IEM Graz for pointing me to the symplectic integration methods, which are essential for audifying ODEs, as they help to ensure numeric stability in the long run (e.g. to avoid drifts of oscillations that are mathematically expected to be regular). See the chapter on integration in his dissertation [2], pp 135-146. You might also want check David Pirròs optimized ODE compiler named Henri [3]. Big credit also to Nathaniel Virgo who brought up the buffering strategy used in Fb1, which is Fb1_ODE's working horse.
WARNING: The usage of this class is – inherently – highly experimental. Be careful with amplitudes, as always with feedback it can become loud! Sudden blowups might result form the mathematical characteristics of the ODE systems or they might come from parameter changes on which ODEs can react extremely sensitive to, they can also stem from numerical accumulation effects. It is highly recommended to take precautionary measures, e.g. by limiting/distorting operators (tanh, clip, softclip, distort) with the compose option (See Fb1_ODE help Ex.5) and/or external limiting and/or using MasterFX from the JITLibExtensions quark.
NOTE: Fb1_ODE in its plain form (without tMul modulation, use of compose etc.) produces audio data as a numerical integration of an ODE initial value problem, defined by Fb1_ODEdef and a numerical procedure, defined by Fb1_ODEintdef. This of course supposes well-defined ODE systems and it should be kept in mind that the Fb1_ODE framework doesn't perform any mathematical checks regarding the principal existence and uniqueness of a solution of a given Fb1_ODEdef.
NOTE: The convenience of direct definition of the ODE relation comes with the price of a large number of UGens involved. You might want to allow a higher number of UGens with the server option numWireBufs. For a nice workflow I'd recommended to take reduced blockSizes (e.g. 1, 2, 4, 8, 16) while experimenting as compile time is shorter, but once you have finished the design of a SynthDef it might pay going back to blocksize 32 or 64 for runtime efficiency, especially if many kr UGens are involved.
See also: Fb1_ODE, Fb1_ODEintdef, Fb1, Fb1_MSD, Fb1_SD, Fb1_Lorenz, Fb1_Hopf, Fb1_HopfA, Fb1_HopfAFDC, Fb1_VanDerPol, Fb1_Duffing
References:
[1] Trefethen, Lloyd N.; Birkisson Ásgeir; Driscoll, Tobin A. (2017):
Exploring ODEs. SIAM - Society for Industrial and Applied Mathematics.
Free download from: https://people.maths.ox.ac.uk/trefethen/Exploring.pdf
[2] Pirrò, David (2017). Composing Interactions. Dissertation.
Institute of Electronic Music and Acoustics, University of Music and Performing Arts Graz.
Free download from: https://pirro.mur.at/media/pirro_david_composing_interactions_print.pdf
[3] https://git.iem.at/davidpirro/henri
Creation / Class Methods
*new (name, function, t0, y0, outScale = 1, diffOutScale = 1)
Creates a new Fb1_ODEdef object and stores it in a Dictionary for further usage with Fb1_ODE.
name - The name of the ODE, expects a Symbol.
Default Fb1_ODEdefs cannot be overwritten.
function - The implementation of the function F which describes the ODE given as
Y'(t) = F(t, Y(t)) where Y, F can be a single-valued or vector-valued.
It must take time as first and a system state (possibly an array) as second parameter
and optional further args that can also be arrays.
The system state is expected to have the size of the system (given by y0's size),
the function's return value must also equal this size.
If the system size is greater than 1 the components of the output array
should be rather written as Functions, as this optimizes the compile process
with symplectic integration procedures, it isn't compulsory though.
See examples below and in the source file Fb1_ODEdef.sc.
t0 - Number. Default initial time value.
y0 - Default initial value of the ODE.
Expects number or array and determines the size of the system, function's second arg
and return value must be of same size.
outScale - Number that determines the default multiplier for the integration signal
produced by Fb1_ODE when the latter's withOutScale flag is set to true (default).
Defaults to 1. Especially thought for systems like Lorenz which produce high levels
with standard parameters, then it makes sense to set outScale to a value smaller than 1.
WARNING: outScale / diffOutScale / withOutScale do not implement a
general safety net functionality, withOutScale's default value true does by no means say
that output is limited. Scaling values of Fb1_ODEdefs can only be average assumptions and can,
under different circumstances, always lead to high out levels.
diffOutScale - Number that determines the default multiplier for the differential signal
produced by Fb1_ODE when the latter's withOutScale flag is set to true (default).
Defaults to 1. Especially thought for systems like Lorenz which produce high levels
with standard parameters, then it makes sense to set diffOutScale to a value smaller than 1.
WARNING: outScale / diffOutScale / withOutScale do not implement a
general safety net functionality, withOutScale's default value true does by no means say
that output is limited. Scaling values of Fb1_ODEdefs can only be average assumptions and can,
under different circumstances, always lead to high out levels.
*at (key)
Returns the Fb1_ODEdef instance of the Symbol key if it exists.
*keys
Returns an array of all keys of currently stored Fb1_ODEdefs.
*postAll
Posts all keys of currently stored Fb1_ODEdefs.
*remove (key)
Removes the Fb1_ODEdef of the Symbol key from the Dictionary.
*reset
Removes all Fb1_ODEdefs other than the predefined ones from the Dictionary.
Instance Methods
next (intType, t, y, dt, args)
Method for language-side integration, gives next value based on previous data.
intType - Integration type.
Expects one of the Symbols, for which procedures are stored with Fb1_ODEintdef.
The use of symplectic procedures (with prefix 'sym', like 'sym2') is highly recommended.
For more accurate integration you can try symplectic procedures of
higher order like \sym4, \sym8 etc. Multi-step procedures are not implemented, remaining:
Symplectic:
\sym2, \sym2_d, \sym4, \sym4_d, \sym6, \sym6_d, \sym8, \sym8_d,
\sym12, \sym12_d, \sym16, \sym16_d, \sym32, \sym32_d, \sym64, \sym64_d
Euler:
\eu, \eu_d, \eum, \eum_d, \eui, \eui_d
Prediction-Evaluation-Correction:
\pec, \pece, \pecec, \pecece
Runge-Kutta:
\rk3, \rk3_d, \rk3h, \rk3h_d, \rk4, \rk4_d
t - Previous time, expects Number.
y - Previous state, expects Number or Array of system size.
dt - Number. Time delta.
args - List of additional args to be passed to Fb1_ODEdef's function.
nextN (n, intType, t, y, dt, args, withTime = false, includeStart = true)
Method for language-side integration, gives an array of next n values (arrays) based on previous data.
Last values (arrays) are stored at last position. For intTypes with sizeFactor 2 the differential is included.
See Examples.
n - Integer. Number of next values (resp. value arrays) to be calculated.
intType - Integration type.
Expects one of the Symbols, for which procedures are stored with Fb1_ODEintdef.
The use of symplectic procedures (with prefix 'sym', like 'sym2') is highly recommended.
For more accurate integration you can try symplectic procedures of
higher order like \sym4, \sym8 etc. Multi-step procedures are not implemented, remaining:
Symplectic:
\sym2, \sym2_d, \sym4, \sym4_d, \sym6, \sym6_d, \sym8, \sym8_d,
\sym12, \sym12_d, \sym16, \sym16_d, \sym32, \sym32_d, \sym64, \sym64_d
Euler:
\eu, \eu_d, \eum, \eum_d, \eui, \eui_d
Prediction-Evaluation-Correction:
\pec, \pece, \pecec, \pecece
Runge-Kutta:
\rk3, \rk3_d, \rk3h, \rk3h_d, \rk4, \rk4_d
t - Previous time, expects Number.
y - Previous state, expects Number or Array of system size.
dt - Number. Time delta.
args - List of additional args to be passed to Fb1_ODEdef's function.
withTime - Boolean. Determines if integrated time should be included.
Defaults to false.
includeStart - Boolean. Determines if start value(s) should be included.
Defaults to true.
name
Getter for the Fb1_ODEdef's name.
function
Getter for the Fb1_ODEdef's function.
t0
Getter for the Fb1_ODEdef's t0.
y0
Getter for the Fb1_ODEdef's y0.
outScale
Getter for the Fb1_ODEdef's outScale.
diffOutScale
Getter for the Fb1_ODEdef's diffOutScale.
size
Getter for the Fb1_ODEdef's system size.
Examples 1) Defining new ODEs
// reboot with reduced blockSize
(
s = Server.local;
Server.default = s;
s.options.blockSize = 8;
s.reboot;
)
Ex.1a) Extending the harmonic oscillator
// This seems to be an interesting strategy.
// We can e.g. try to multiply one of its equations with a term near 1,
// so start with rather small k and investigate
// y'(t) = w(t)
// w'(t) = -y(t) * (1 + (k * w(t))
(
Fb1_ODEdef(\harmonic_ext_1, { |t, y, k|
[
y[1],
y[0].neg * (1 + (k * y[1]))
]
}, 0, [0, 1], 1, 1);
)
// brassy sound
(
x = {
var sig = Fb1_ODE.ar(\harmonic_ext_1,
[1], 1500, 0, [0, 1],
) * 0.1;
sig
}.play
)
x.release
// with oscillating k the system tends to become unstable
// but with softclip per sample it can be kept
// (see chapter 'compose' option in Fb1_ODE help)
(
x = {
var sig = Fb1_ODE.ar(\harmonic_ext_1,
// k oscillates between 1 and 3
[SinOsc.ar(50).lincurve(-1, 1, 1, 3, 2)],
1500, 0, [0, 1],
compose: \softclip
) * 0.1;
sig
}.play
)
x.release
// an oscillation between 1 and higher values totally changes the spectrum
(
x = {
var sig = Fb1_ODE.ar(\harmonic_ext_1,
// k oscillates between 1 and 10
[SinOsc.ar(50).lincurve(-1, 1, 1, 10, 2)],
1500, 0, [0, 1],
compose: \softclip
) * 0.1;
sig
}.play
)
x.release
// play with the two states
(
x = {
var sig = Fb1_ODE.ar(\harmonic_ext_1,
// upper oscillation bound for k oscillates itself between 2 and 15
[SinOsc.ar(50).lincurve(-1, 1, 1, SinOsc.ar(SinOsc.ar(0.2).exprange(0.2, 15)).range(2, 15), 2)],
1500, 0, [0, 1],
compose: \softclip
) * 0.1;
sig
}.play
)
x.release
Ex.1b) Extending exponential decay
// exponential decay is described by the equation
// y'(t) = -y(t)
// an oscillating decay can e.g. be got by
// y'(t) = -y(t) * sin(t)
// the analytic solution includes a log of the sine,
// so we get more partials
(
Fb1_ODEdef(\exp_decay_raw, { |t, y|
y.neg * sin(t)
}, 0, 1, 1, 1);
)
(
x = {
var sig = Fb1_ODE.ar(\exp_decay_raw,
tMul: 100 * 2pi,
compose: \softclip
) ! 2;
Line.kr(dur: 10, doneAction: 2);
sig
}.play
)
x.release
// multiplication with a second sine with multiplied time leads to strange and interesting results
(
Fb1_ODEdef(\exp_decay_extended, { |t, y, k|
y.neg * (sin(t) * sin(k * t))
}, 0, 1, 1, 1);
)
// ATTENTION: danger of blowup, can be reduced with softclip composition per sample
// constant values lead to ring modulation-like effects ...
(
x = {
var sig = Fb1_ODE.ar(\exp_decay_extended,
[2.7], 100 * 2pi, 0, 1,
compose: \softclip
);
Line.kr(dur: 10, doneAction: 2);
sig ! 2
}.play
)
x.release
// ... whereas modulations produce more complex changing spectra
(
x = {
var sig = Fb1_ODE.ar(\exp_decay_extended,
[SinOsc.ar(120).range(3, 3.01)], 100 * 2pi, 0, 1,
compose: \softclip
);
Line.kr(dur: 10, doneAction: 2);
sig ! 2
}.play
)
x.release
// for decorrelated stereo we can expand to two independent equations
// k should be of size 2
(
Fb1_ODEdef(\exp_decay_extended_2, { |t, y, k|
[
y[0].neg * (sin(t) * sin(k[0] * t)),
y[1].neg * (sin(t) * sin(k[1] * t))
]
}, 0, [1, 1], 1, 1);
)
(
x = {
var sig = Fb1_ODE.ar(\exp_decay_extended_2,
[SinOsc.ar(120).range([3, 3.01], [3.01, 3.02])], 100 * 2pi, 0, [1, 1],
compose: \softclip
);
Line.kr(dur: 10, doneAction: 2);
sig
}.play
)
x.release
Ex.2) Language-side integration
// Possible though not optimized.
// For longer sections it's probably better to employ Fb1 ODE solvers,
// store audio in a buffer and load it into a float array,
// for quick tests it might be useful though.
// integrate begin of a Lorenz system, last array is last state, so flop for plot
(
a = Fb1_ODEdef.at(\Lorenz).nextN(
n: 1000,
intType: \sym2,
t: 0,
y: [1, 1, 1],
dt: 1/100,
args: [30, 12, 2]
);
a.flop.plot
)