Fb1_ODE | SuperCollider 3.13.0 Help

Description

Pseudo ugen for integrating / audifying ordinary (systems of) differential equations with initial values in realtime, based on the Fb1 single sample feedback class. It makes use of Fb1_ODEdef and Fb1_ODEintdef, containers for ODEs and numerical solution methods, which provide an interface for adding new ODE systems or integration methods interactively. There exists a huge number of ODE models in different areas, e.g. have a look at the SIAM publication Exploring ODEs [1], it approaches the topic from an experimental point of view, which is quite related to the demands of practical audification / synthesis / processing. And of course it's fun to define your own ODEs, see this help file and Fb1_ODEdef.

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: Especially with self-defined ODEs 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 Examples 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.

References

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
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
https://git.iem.at/davidpirro/henri

Class Methods

Fb1_ODE.new(name, argList, tMul: 1, t0, y0, intType: 'sym2', compose, composeArIn, dt0, argList0, init_intType: 'sym8', withOutScale: true, withDiffChannels: false, withTimeChannel: false, blockSize, graphOrderType: 1, leakDC: true, leakCoef: 0.995)

Creates a new Fb1_ODE ar object.

Arguments:

name	The name of the ODE defined and stored via Fb1_ODEdef. Can be one of the predefined ODEs, for which wrappers exist, or a self-defined one. Expects a Symbol.
argList	The argList to be passed to the ODE. Can contain ar or kr UGens and SequenceableCollections thereof.
tMul	Time multiplier, which determines velocity of proceeding in the dynamic system. The default value 1 means that the step delta of time equals 1 / sample duration. For getting audible oscillations you might, depending on the ODE definition, have to pass much higher values. Can also be a kr / ar UGen and might also be negative.
t0	Initial time. Expects Number. If no value is passed, the default value of the referred Fb1_ODEdef is taken, usually 0.
y0	Initial value of the ODE. Expects number or array of correct size, which is determined by the Fb1_ODEdef. If no value is passed, the default value of the referred Fb1_ODEdef is taken.
intType	Integration type. Expects one of the Symbols, for which procedures are stored with Fb1_ODEintdef. The use of symplectic procedures (with prefix 'sym') is highly recommended. Defaults to \sym2. For more accurate integration you can try symplectic procedures of higher order like \sym4, \sym8 etc. Families of integration procedures: 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 Adams-Bashforth: \ab2, \ab3, \ab4, \ab5, \ab6 Adams-Bashforth-Moulton: \abm21, \abm22, \abm32, \abm33, \abm43, \abm44, \abm54, \abm55, \abm65, \abm66
compose	Operator(s) / Function(s) to be applied to the system value on a per-sample base. This of course blurs the numeric procedure but can be used for containing or in a creative way. Can be an operator Symbol, a Function or an arbitrarily mixed SequenceableCollection thereof. The Functions are in any case expected to take an array (the system state) as first argument. If only one Function is given it must output an array of same (system) size. Within a SequenceableCollection a Function must output a value of size 0. Within a SequenceableCollection an operator Symbol is applied only to the corresponding component. A Function can optionally take a second argument which is for ar UGens passed via composeArIn.
composeArIn	ar UGen or SequenceableCollection thereof or a SequenceableCollection that can contain both. This is the way to use ar UGens within a Function passed to compose. UGens are passed to the Function's second argument.
dt0	First time delta in seconds to be used for the language-side calculation of first values of a multi-step intType procedure. This will mostly be irrelevant as (single-step) symplectic procedures are to be preferred. In case of a multi-step procedure a dt0 value will be derived from the default server's properties (sample duration * tMul).
argList0	Initial argList value(s) to be used for the language-side calculation of first values of a multi-step intType procedure. This will mostly be irrelevant as (single-step) symplectic procedures are to be preferred. If no UGens are passed to argList, values will be assumed from passed argList Numbers.
init_intType	Integration type for language-side calculation of first values of a multi-step intType procedure. This will mostly be irrelevant as (single-step) symplectic procedures are to be preferred. Defaults to \sym8.
withOutScale	Boolean. Determines if the Fb1_ODEdef's default scaling values for integration and differential signals should be applied to the output. Defaults to true. WARNING: withOutScale does not implement a general safety net functionality, withOutScale's default value true does by no means say that output is limited. The default scaling values of predefined Fb1_ODEdefs are average assumptions and can, under different circumstances, always lead to high out levels.
withDiffChannels	Boolean. Determines if channel(s) with differential value(s) should be returned. This is only applicable with integration types with a sizeFactor == 2, which means that for every sample the values of the differential are buffered. As by default this is not the case for all predefined numeric procedures, there exist variants with appendix '_d', e.g. \sym2_d. Defaults to false.
withTimeChannel	Boolean. Determines if accumulated time is output in an additional channel. WARNING: with constant tMul it produces an ascending DC which is not affected by leakDC if this is set to true. Might especially be of interest if time is modulated. Defaults to false.
blockSize	Integer, this should be the server blockSize. If no Integer is passed, the current default server's blockSize is assumed (in contrast to Fb1). So explicitely passing a blockSize should only be necessary in special cases, e.g. when compiling before booting. However for a pleasant workflow for ar usage it's recommended to use Fb1_ODE with a reduced server blockSize in order to reduce SynthDef compile time.
graphOrderType	0, 1 or 2. Determines if topological order of generated BufRd and BufWr instances in the SynthDef graph is forced by additional UGens. Type 0: forced graph order is turned off. Type 1 (default): graph order is forced by summation and <!. Type 2: graph order is forced by <! operators only. Default 1 is recommended, but with CPU-intense SynthDefs it might be worth trying it with the value 0. This saves a lot of UGens though there might be exceptional cases with different results. Type 2 can shorten the SynthDef compilation time for certain graphs with a large number of UGens, which can be lengthy with type 1. However, CPU usage doesn't directly correspond to the number of UGens.
leakDC	Boolean. Determines if a LeakDC is applied to the output (except time channel). Defaults to true.
leakCoef	Number, the leakDC coefficient. Defaults to 0.995.

Fb1_ODE.ar(name, argList, tMul: 1, t0, y0, intType: 'sym2', compose, composeArIn, dt0, argList0, init_intType: 'sym8', withOutScale: true, withDiffChannels: false, withTimeChannel: false, blockSize, graphOrderType: 1, leakDC: true, leakCoef: 0.995)

Equivalent to *new.

Fb1_ODE.kr(name, argList, tMul: 1, t0, y0, intType: 'sym4', compose, dt0, argList0, init_intType: 'sym8', withOutScale: true, withDiffChannels: false, withTimeChannel: false, graphOrderType: 1, leakDC: true, leakCoef: 0.995)

Creates a new Fb1_ODE kr object.

Arguments:

name	The name of the ODE defined and stored via Fb1_ODEdef. Can be one of the predefined ODEs, for which wrappers exist, or a self-defined one. Expects a Symbol.
argList	The argList to be passed to the ODE. Can contain kr UGens and SequenceableCollections thereof.
tMul	Time multiplier, which determines velocity of proceeding in the dynamic system. The default value 1 means that the step delta of time equals 1 / control duration. For getting audible oscillations you might, depending on the ODE definition, have to pass much higher values. Can also be a kr UGen and might also be negative.
t0	Initial time. Expects Number. If no value is passed, the default value of the referred Fb1_ODEdef is taken, usually 0.
y0	Initial value of the ODE. Expects number or array of correct size, which is determined by the Fb1_ODEdef. If no value is passed, the default value of the referred Fb1_ODEdef is taken.
intType	Integration type. Expects one of the Symbols, for which procedures are stored with Fb1_ODEintdef. The use of symplectic procedures (with prefix 'sym') is highly recommended. Defaults to \sym4. For more accurate integration you can try symplectic procedures of higher order like \sym4, \sym8 etc. Families of integration procedures: 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 Adams-Bashforth: \ab2, \ab3, \ab4, \ab5, \ab6 Adams-Bashforth-Moulton: \abm21, \abm22, \abm32, \abm33, \abm43, \abm44, \abm54, \abm55, \abm65, \abm66
compose	Operator(s) / Function(s) to be applied to the system value on a per-sample base. This of course blurs the numeric procedure but can be used for containing or in a creative way. Can be an operator Symbol, a Function or an arbitrarily mixed SequenceableCollection thereof. The Functions are in any case expected to take an array (the system state) as first argument. If only one Function is given it must output an array of same (system) size. Within a SequenceableCollection a Function must output a value of size 0. Within a SequenceableCollection an operator Symbol is applied only to the corresponding component.
dt0	First time delta in seconds to be used for the language-side calculation of first values of a multi-step intType procedure. This will mostly be irrelevant as (single-step) symplectic procedures are to be preferred. In case of a multi-step procedure a dt0 value will be derived from the default server's properties (control duration * tMul).
argList0	Initial argList value(s) to be used for the language-side calculation of first values of a multi-step intType procedure. This will mostly be irrelevant as (single-step) symplectic procedures are to be preferred. If no UGens are passed to argList, values will be assumed from passed argList Numbers.
init_intType	Integration type for language-side calculation of first values of a multi-step intType procedure. This will mostly be irrelevant as (single-step) symplectic procedures are to be preferred. Defaults to \sym8.
withOutScale	Boolean. Determines if the Fb1_ODEdef's default scaling values for integration and differential signals should be applied to the output. Defaults to true. WARNING: withOutScale does not implement a general safety net functionality, withOutScale's default value true does by no means say that output is limited. The default scaling values of predefined Fb1_ODEdefs are average assumptions and can, under different circumstances, always lead to high out levels.
withDiffChannels	Boolean. Determines if channel(s) with differential value(s) should be returned. This is only applicable with integration types with a sizeFactor == 2, which means that for every sample the values of the differential are buffered. As by default this is not the case for all predefined numeric procedures, there exist variants with appendix '_d', e.g. \sym2_d. Defaults to false.
withTimeChannel	Boolean. Determines if accumulated time is output in an additional channel. WARNING: with constant tMul it produces an ascending DC which is not affected by leakDC if this is set to true. Might especially be of interest if time is modulated. Defaults to false.
graphOrderType	0, 1 or 2. Determines if topological order of generated BufRd and BufWr instances in the SynthDef graph is forced by additional UGens. Type 0: forced graph order is turned off. Type 1 (default): graph order is forced by summation and <!. Type 2: graph order is forced by <! operators only. Default 1 is recommended, but with CPU-intense SynthDefs it might be worth trying it with the value 0. This saves a lot of UGens though there might be exceptional cases with different results. Type 2 can shorten the SynthDef compilation time for certain graphs with a large number of UGens, which can be lengthy with type 1. However, CPU usage doesn't directly correspond to the number of UGens.
leakDC	Boolean. Determines if a LeakDC is applied to the output (except time channel). Defaults to true.
leakCoef	Number, the leakDC coefficient. Defaults to 0.995.

Inherited class methods

Instance Methods

Inherited instance methods

Examples 1: Proof of concept

Ex. 1a: The harmonic oscillator

// A sine/cos function is the solution of the second order differential equation

// y''(t) = -y(t)

// with the standard substitution ...

// w(t) = y'(t)

// ... this can be reformulated as a system of two equations of first order

// y'(t) = w(t)
// w'(t) = -y(t)

// In general an ODE can be written like this
// where Y and F are meant to be a vector-valued functions,
// F given and Y to be found

// Y'(t) = F(t, Y)

// For the corresponding Fb1_ODEdef the Function (F) must take time and - as here we have size 2 - 
// an array as arguments plus optional further arguments and 
// output an array of same size, which can be implemented by the evaluation of:

(
Fb1_ODEdef(\harmonic, { |t, y|
    [y[1], y[0].neg]
}, 0, [0, 1], 1, 1); // default init values t0, y0 and scaling factors for output
)

// The stereo output is a pair of 90-degree-shifted sine waves.
// As one wave length would be 2pi seconds we multiply time by a factor 500 * 2pi to get 500 Hz.
// Changed blockSize is detected automatically.

x = { Fb1_ODE.ar(\harmonic, tMul: 500 * 2pi) * 0.1 }.play

x.release

// as with Fb1, method 'new' is equivalent to 'ar'

x = { Fb1_ODE(\harmonic, tMul: 500 * 2pi) * 0.1 }.play

x.release

// a default frequency could equivalently also be built into the model

(
Fb1_ODEdef(\harmonic_2, { |t, y|
    [y[1], y[0].neg] * 500 * 2pi
}, 0, [0, 1], 1, 1); // default init values t0, y0 and scaling factors for output
)

x = { Fb1_ODE.ar(\harmonic_2) * 0.1 }.play

x.release

Ex. 1b: Extending the harmonic oscillator to the mass-spring-damper model

// The mass-spring-damper model with externally applied force
// is described by the second order differential equation

// y''(t) * mass = f(t) - (dampen * y'(t)) - (spring * y(t))

// where dampen means the dampening factor and spring the spring stiffness.
// Again standard substitution

// w(t) = y'(t)

// leads to the system of two first order equations

// y'(t) = w(t)
// w'(t) = (f(t) - (dampen * w(t)) - (spring * y(t))) / mass

// In the source code of Fb1_ODEdef.sc the corresponding Fb1_ODEdef looks like this:
// the Function brackets within the array are not absolutely necessary, 
// but compile process is faster.

// don't evaluate, already stored

Fb1_ODEdef(\MSD, { |t, y, f = 0, mass = 1, spring = 1, dampen = 0|
    [
        { y[1] },
        { f - (dampen * y[1]) - (spring * y[0]) / mass }
    ]
}, 0, [0, 0], 1, 1);

// Here the resulting oscillation, which describes the position, is used for FM,
// the deflection converges to a value of 0.004, so take no LeakDC.

// LR difference results from slightly different linear mapping

(
x = {
    // we need extreme values for fast oscillation
    var f = 0.2, mass = 0.001, spring = 50, dampen = 0.0005,

// for t0, y0 default values are taken
        sig = Fb1_ODE.ar(\MSD, 
            [f, mass, spring, dampen],
            leakDC: false
        );
    SinOsc.ar(sig[0].linlin(0, 0.005, [100, 101], 700), 0, 0.1)
}.play
)

x.release

// for audible sound production take a time multiplier and a scaling factor
// and get a decaying sine, leakDC now true by default

(
x = {
    var f = 0.2, mass = 0.001, spring = 50, dampen = 0.0005,
        sig = Fb1_ODE.ar(\MSD, [f, mass, spring, dampen], tMul: 10);
    Line.kr(dur: 10, doneAction: 2);
    sig[0] * 20
}.play
)

x.release

// as MSD is predefined together with the wrapper class Fb1_MSD
// one can equivalently write

(
x = {
    var f = 0.2, mass = 0.001, spring = 50, dampen = 0.0005,
        sig = Fb1_MSD.ar(f, mass, spring, dampen, tMul: 10);
    Line.kr(dur: 10, doneAction: 2);
    sig[0] * 20
}.play
)

x.release

Ex. 2: Numeric integration type

Examples 3: Modulations

All optional arguments of the Fb1_ODE can take UGens at audio or control rate, tMul can be modulated at both rates as well.

Ex. 3a: Argument modulation

// additional oscillation of mass, this is already a quite extreme blurring
// operations of such kind can easily result in derailing the oscillaton

// here it's interesting to hear quite a difference between ar and kr variants

// ar modulation

(
x = {
    var f = 0.2, mass = 0.001, spring = 50, dampen = 0.0005,
        sig = Fb1_ODE.ar(\MSD, 
            [f, mass * LFTri.ar(200).range(0.1, 2), spring, dampen],
            leakDC: false
        );
    SinOsc.ar(sig[0].linlin(0, 0.005, [100, 101], 700), 0, 0.1)
}.play
)
    
x.release

// much less smooth with kr

(
x = {
    var f = 0.2, mass = 0.001, spring = 50, dampen = 0.0005,
        sig = Fb1_ODE.ar(\MSD, 
            [f, mass * LFTri.kr(200).range(0.1, 2), spring, dampen],
            leakDC: false
        );
    SinOsc.ar(sig[0].linlin(0, 0.005, [100, 101], 700), 0, 0.1)
}.play
)

x.release

// due to Fb1's design, above kr modulator is not interpolated
// an interpolated kr variant though gives an even more blurred result

(
x = {
    var f = 0.2, mass = 0.001, spring = 50, dampen = 0.0005,
        sig = Fb1_ODE.ar(\MSD, 
            [f, mass * K2A.ar(LFTri.kr(200).range(0.1, 2)), spring, dampen],
            leakDC: false
        );
    SinOsc.ar(sig[0].linlin(0, 0.005, [100, 101], 700), 0, 0.1)
}.play
)

x.release

Ex. 3b: Time modulation

Ex. 4: Initial values

Examples 5: The 'compose' argument

This allows to build an additional Function into the Fb1 feedback loop. It is applied to every array of samples that is result and next input of the numeric integration procedure. Obviously then correct integration of the ODE cannot be performed anymore, but the option has still its value.

Ex. 5a: Handling instable systems

// mass spring damper with extra term y[0] * y[1] -
// out of interest, no physical argument for that

(
Fb1_ODEdef(\MSD_2, { |t, y, f = 0, mass = 1, spring = 1, dampen = 0|
    [
        y[1], 
        f - (dampen * y[1]) - (spring * y[0]) + (y[0] * y[1]) / mass        
    ]
}, 0, [0, 0], 1, 1);
)

// System derails after some seconds and produces nans.

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var sig = Fb1_ODE.ar(\MSD_2, [f, mass, spring, dampen], leakDC: false).poll;
    LeakDC.ar(SinOsc.ar(sig[0].linlin(0, 0.006, [50, 50.1], 700), 0, 0.1))
}.play
)

x.release

// Keeping the system alive with clipping at LFO-controlled frequency.
// Note that the compose Function is passed an array as size of MSD_2 is 2,
// clip2 applied to an array again returns an array which is necessary,
// the Function must preserve the system size.

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var lfo = LFDNoise3.kr(2).exprange(0.2, 1500);
    var sig = Fb1_ODE.ar(\MSD_2, 
        [f, mass, spring, dampen], 
        leakDC: false,
        compose: { |y| y.clip2(lfo) }        
    );
    SinOsc.ar(sig[0].linlin(0, 0.006, [50, 50.1], 700).poll, 0, 0.1)
}.play
)

x.release

// if an operator is passed via a Symbol, it applies to all channels

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var sig = Fb1_ODE.ar(\MSD_2, 
        [f, mass, spring, dampen], 
        leakDC: false,
        compose: \softclip        
    );
    SinOsc.ar(sig[0].linlin(0, 0.006, [50, 50.1], 700).poll, 0, 0.1)
}.play
)

x.release

Ex. 5b: Other options of the 'compose' argument

// You might want to pass ar signals to the composition Function.
// This has to be done via the composeArIn argument, 
// its UGens are then passed as second argument to the compose Function.

// This somewhat odd way of passing is necessary as feeding ar signals into Fb1 is not trivial
// and must also be done via an extra 'in' argument - which in turn is used by Fb1_ODE.

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var mod = SinOsc.ar(SinOsc.ar(LFDNoise3.ar(0.1).range(50, 1000)).range(10, 5000), 0, 0.1, 1);
    var sig = Fb1_ODE.ar(\MSD,
        [f, mass, spring, dampen],
        leakDC: false,
        compose: { |y, in| (y * in).tanh },
        composeArIn: mod
    );
    SinOsc.ar(sig[0].linlin(0, 0.006, [50, 50.1], 700), 0, 0.07)
}.play
)

x.release

// note that above operation is totally different from applying the ring modulation + tanh
// after the integration, here the feedback process is much simpler:

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var mod = SinOsc.ar(SinOsc.ar(LFDNoise3.ar(0.1).range(50, 1000)).range(10, 5000), 0, 0.1, 1);
    var sig = (Fb1_ODE.ar(\MSD,
        [f, mass, spring, dampen],
        leakDC: false
    ) * mod).tanh;
    SinOsc.ar(sig[0].linlin(0, 0.006, [50, 50.1], 700), 0, 0.07)
}.play
)

x.release

// The compose arg can also take an array of Functions and/or operators,
// the composeArIn arg can take an array of UGens,
// the compose Function(s) can refer to the components of the passed arrays.

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var mod = SinOsc.ar(SinOsc.ar(LFDNoise3.ar(0.1 ! 2).range(50, 100)).range(10, 5000), 0, 0.1, 1);
    var sig = Fb1_ODE.ar(\MSD,
        [f, mass, spring, dampen],
        leakDC: false,
        compose: [
            { |y, in| (y[0] * in[0]).softclip },
            { |y, in| ((y[0] + y[1]) * in[1]).softclip }
        ],
        composeArIn: mod // mod is a stereo signal !
    );
    SinOsc.ar(sig[0..1].linlin(0, 0.006, 50, LFDNoise1.ar(0.1).range(700, 2000)), 0, 0.05)
}.play
)

x.release

// composeArIn's array can contain arrays also

(
x = {
    var f = 0.5, mass = 0.001, spring = 150, dampen = 0.001;
    var mod0 = Saw.ar(LFDNoise1.ar(1).exprange(3, 3000), 0.05);
    var mod1 = SinOsc.ar(SinOsc.ar(LFDNoise3.ar(1 ! 2).range(1, 20)).range(1, 3000), 0, 0.1, 1);
    var sig = Fb1_ODE.ar(\MSD,
        [f, mass, spring, dampen],
        leakDC: false,
        compose: [
            { |y, in| (y[0] * in[1][0]).softclip },
            { |y, in| ((y[0] + y[1]) * in[1][1] + in[0]).softclip }
        ],
        composeArIn: [mod0, mod1] // mod1 is itself a stereo signal !
    );
    SinOsc.ar(sig[0..1].linlin(0, 0.006, 50, LFDNoise1.ar(0.1).range(700, 1500)), 0, 0.05);
}.play
)

x.release

Examples 6: Channels with additional information

Ex. 6a: The differential channels

// The optional differential channel(s) contain the differential(s) at time t,
// in other terms: the value(s) F(t, Y) of the system
// Y'(t) = F(t, Y)

// For many numeric integration methods these values are buffered anyway,
// but not for the symplectic procedures.
// Therefore there exist variants with suffix '_d' which can be used in that case

// an oscillation produced by a Hopf ODE (see also Fb1_Hopf)
    
(
x = {
    Fb1_ODE.ar(\Hopf,
        [1, 2, 1],
        1000,
        intType: \sym2_d, // no difference to sym2 in the first two channels
    ) * 0.5
}.play
)

x.release
    
// The system has size two, so with differential we get 4 out channels,
// take the differential only, it gives a different sound colour (more partials).
    
(
x = {
    Fb1_ODE.ar(\Hopf,
        [1, 2, 1],
        1000,
        intType: \sym2_d,
        withDiffChannels: true,
    )[2..3] * 0.5
}.play
)

x.release
    
    
// The differential channels are basically buffering the slope of the 
// original (integration/solution) signal as a by-product, 
// slope could also be calculated with the Slope UGen.
// Note that we scaled time by a factor 1000, so the result of Slope
// has to be divided by the same factor.
    
(
x = {
    Slope.ar(
        Fb1_ODE.ar(\Hopf,
            [1, 2, 1],
            1000,
            intType: \sym2_d,
            withDiffChannels: true
        )[0..1] * 0.5
    ) / 1000
}.play
)

x.release

Ex. 6b: The time channel

// tMul also multiplies the integrated time, so here we proceed with 500 * 2pi per second ...

(
Fb1_ODEdef(\harmonic, { |t, y|
    [y[1], y[0].neg]
}, 0, [0, 1], 1, 1); 
)
    
(
x = { 
    var sig = Fb1_ODE.ar(\harmonic, 
        tMul: 500 * 2pi,
        withTimeChannel: true    
    );
    sig[2].poll;
    sig[0..1] * 0.1 // don't want to have time as audio output !
}.play
)
    
x.release
    
        
// ... whereas here time is scaled by the ODE    
    
(
Fb1_ODEdef(\harmonic_2, { |t, y|
    [y[1], y[0].neg] * 500 * 2pi
}, 0, [0, 1], 1, 1); 
)
    
(
x = { 
    var sig = Fb1_ODE.ar(\harmonic_2, 
        tMul: 1,
        withTimeChannel: true    
    );
    sig[2].poll;
    sig[0..1] * 0.1 // don't want to have time as audio output !
}.play
)
    
x.release

// time integration might be of interest with modulations of tMul
    
(
x = {
    var sig = Fb1_ODE.ar(\harmonic_2,
        tMul: SinOsc.kr(0.2).exprange(0.2, 2),
        withTimeChannel: true
    );
    sig[2].poll;
    sig[0..1] * 0.07 // don't want to have time as audio output !
}.play
)

x.release

Ex. 7: The 'withOutScale' argument

// It determines if default scaling parameters for the ODE integration signal and differential(s) should be applied
// In many cases this doesn't make a difference when the default scaling values equal 1.

// Certain ODEs though produce a very high amplitude level with usual standard params, like Lorenz.
// So it makes sense to scale the output down by default.

// don't evaluate, already stored

Fb1_ODEdef(\Lorenz, { |t, y, s = 10, r = 30, b = 2|
    [
        { (y[1] - y[0]) * s },
        { (r - y[2]) * y[0] - y[1] },
        { (y[0] * y[1]) - (b * y[2]) }
    ]
}, 0, [1, 1, 1], 0.01, 0.01); // scaling parameters 0.01

// So there's no problem to use Lorenz with standard params as
// 'withOutScale' defaults to true

(
x = {
    var s = 10, r = 30, b = 3, tMul = MouseX.kr(50, 150, 1).poll,
        sig;
    sig = Fb1_ODE.ar(\Lorenz, [s, r, b], tMul);
    sig[0..1]
}.play;
)

x.release

// WARNING: it's not recommended to set 'withOutScale' to false
// as in cases like this you get very high levels !
// So here we have to reduce the level outside the Fb1_ODE

(
x = {
    var s = 10, r = 30, b = 3, tMul = MouseX.kr(50, 150, 1).poll,
        sig;
    sig = Fb1_ODE.ar(\Lorenz, [s, r, b], tMul, withOutScale: false);
    sig[0..1] / 100
}.play;
)

x.release

// WARNING: 'withOutScale' does not implement a general safety net functionality, 
// withOutScale's default value true does by no means say that output is limited. 
// The default scaling values of predefined Fb1_ODEdefs are average assumptions and can, 
// under different circumstances, always lead to high out levels.

Examples 8: ODEs from mechanics

Ex. 8a: The driven pendulum

Ex. 8b: The reduced two body problem

// Two planets are moving around the common barycenter, thus the system can be simplified
// to an ODE which describes the changes of the the relative vector between the masses.

// E.g. see https://evgenii.com/blog/two-body-problem-simulator

(
Fb1_ODEdef(\TwoBodyReduced, { |t, y, m1, m2|
    var v = -1 - (m1 / m2);
    var r = y[0..1].squared.sum.sqrt;
    [
        { y[2] },
        { y[3] },
        { v * y[0] / r.cubed },
        { v * y[1] / r.cubed }
    ]
}, 0, [0, 0, 0, 0], 1, 1);
)

// It turns out that this system is numerically very sensitive.
// Standard procedures are failing soon: here classical Runge-Kutta 4rth order
// leads to a glissando before it collapses.
// Even 2nd order symplectic shows movement of wave forms, 
// though it keeps frequency and preserves volume

// here we take only the first two components which represent the plane coordinates
// of the difference vector

// sym8 still shows a slow morphing of wave forms
// you can try with \sym64, very inefficient, but even more steady

(
x = {
    var m = [1, 1.5];
    Fb1_ODE.ar(\TwoBodyReduced, m,
        100, 0, [0, 0.2, 2, 1],
        intType: \sym8
    )[0..1]
}.play
)

s.scope

x.release

Examples 9: ODEs from population dynamics

Ex. 9a: The Lotka-Volterra equations

Ex. 9b: The Hastings-Powell equations

// This is an extension of the Lotka-Volterra model and
// considers the fact that rabbits must eat too.
// It is described in the article

// Alan Hastings and Thomas Powell (1991). 
// "Chaos in a Three-Species Food Chain" Ecology 72(3): 896-903

// and on page 165 of the book "Exploring ODEs"

// Trefethen, Lloyd N.; Birkisson, Ásgeir; Driscoll, Tobin A. (2017):
// Exploring ODEs. SIAM - Society for Industrial and Applied Mathematics
// Free download from:
// http://people.maths.ox.ac.uk/trefethen/Exploring.pdf

// We take over initial values and parameters from the example given there:

(
Fb1_ODEdef(\CarrotsRabbitsFoxes, { |t, y, a1 = 5, a2 = 0.1,
    b1 = 0.1, b2 = 2, d1 = 0.4, d2 = 0.01|
    var p = a1 * y[0] / (1 + (b1 * y[0])) * y[1];
    var q = a2 * y[1] / (1 + (b2 * y[1])) * y[2];
    [
        { y[0] * (1 - y[0]) - p }, // carrots
        { p - q - (d1 * y[1])  }, // rabbits
        { q - (d2 * y[2])  } // foxes
    ]
}, 0, [0.4, 1, 9], 1, 1);
)

// In stereo we hear population oscillations of carrots and rabbits,
// it takes some fractions of a second at begin after a state is reached.
// Check with mouse and see post window:
// b1 goes through a zone of chaotic behaviour between 2.5 and 3.5

(
x = {
    var a1 = 5, a2 = 0.1, b1 = MouseX.kr(2.5, 3.5).poll, b2 = 2,
        d1 = 0.4, d2 = 0.01, tMul = 5000, amp = 0.2;
    Fb1_ODE.ar(\CarrotsRabbitsFoxes,
        [a1, a2, b1, b2, d1, d2],
        tMul, 0, [0.4, 1, 9]
    )[0..1] * amp
}.play
)

x.release

helpfile source: /home/stefan/.local/share/SuperCollider/downloaded-quarks/miSCellaneous_lib/HelpSource/Classes/Fb1_ODE.schelp
link::Classes/Fb1_ODE::

Fb1_ODE : UGen : AbstractFunction : Object Extension

Fb1_ODE.new(name, argList, tMul: 1, t0, y0, intType: 'sym2', compose, composeArIn, dt0, argList0, init_intType: 'sym8', withOutScale: true, withDiffChannels: false, withTimeChannel: false, blockSize, graphOrderType: 1, leakDC: true, leakCoef: 0.995)

Arguments:

Fb1_ODE.ar(name, argList, tMul: 1, t0, y0, intType: 'sym2', compose, composeArIn, dt0, argList0, init_intType: 'sym8', withOutScale: true, withDiffChannels: false, withTimeChannel: false, blockSize, graphOrderType: 1, leakDC: true, leakCoef: 0.995)

Fb1_ODE.kr(name, argList, tMul: 1, t0, y0, intType: 'sym4', compose, dt0, argList0, init_intType: 'sym8', withOutScale: true, withDiffChannels: false, withTimeChannel: false, graphOrderType: 1, leakDC: true, leakCoef: 0.995)

Arguments:

Fb1_ODE : UGen : AbstractFunction : Object
Extension