Chaos in Faust

[ProgHead]

Uit het blokschema lezen we voor de gesimuleerde DV af dat:

x(t) = A(B(x(t))
x(t) = -B(x(t))
x(t) = -1645333507 * integrator( integrator( arctan(x(t)) ) )

Nu we weten hoe we de Faust "integrator" in een natuurkundig correcte integraal kunnen uitdrukken vinden we:

[imath] x(t) = - 1645333507 \cdot \mathrm{SR}^2 \cdot \int_0^t \int_0^t \arctan(x(\tau)) \, \mathrm{d} \tau \mathrm{d} \tau [/imath]

[imath] \ddot{x}(t) = - 1645333507 \cdot \mathrm{SR}^2 \cdot \arctan(x(t)) [/imath]

[imath] \ddot{x}(t) + 1645333507 \cdot \mathrm{SR}^2 \cdot \arctan(x(t)) = 0 [/imath]

Laat nu:

[imath] \omega_0 = \sqrt{ 1645333507 } \cdot \mathrm{SR} [/imath]

Dan komt er:

[imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x(t)) = 0 [/imath]

Voor SR = 48000 samples/sec vinden we:

[imath] \omega_0 \approx 1.947 \cdot 10^9 \mathrm{rad}/s [/imath]

Voor kleine signalen mogen we dan een golfje van ca. 300 MHz verwachten. Maar onze sample rate is "slechts" 48000 samples/sec...

 

[ProgHead]

De differentiaalvergelijking [imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x(t)) = 0 [/imath] heeft enkel maar periodieke oplossingen. Zie: Does $ \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x(t)) = 0 $ have only periodic solutions?

Dit betekent dat de riedeltjes die worden gegenereerd door het Faust progje ChaosRiedel dat de genoemde differentiaalvergelijking had moeten simuleren digitale artefacten zijn. De oplossingen van de differentiaalvergelijking zelf corresponderen immers enkel met constante tonen, en de geluiden die uit het progje komen lijken daar in de verste verte niet op.

 

[Disharmonic]

Eigenlijk lijkt me dit heel geschikt voor een hybride systeem, een analoge schakeling voor het genereren van een chaotisch signaal en vervolgens met wat digitale toverkunst daar iets bruikbaars van maken dat goed klinkt.

 

[ProgHead]

Ben nog even bezig met een Jan Boerefluitjes bewijs, dat ik ook zelf nog kan volgen , voor het periodieke karakter van de oplossingen van de DV. Is dat eenmaal klaar, dan ga ik uitzoeken hoe we het digitale riedeltjes artefact ook als gewenste feature kunnen gebruiken en aansturen.

 

[ProgHead]

Het bewijs is niet eenvoudig en daarom plaats ik dat hier in meerdere afleveringen.

Stel je een massa m op positie [imath] X(t) = \mathrm{X}_0 \cdot x(t) [/imath] voor aan een veer die een kracht [imath] \mathrm{F}(x) = - \mathrm{F}_0 \cdot \arctan(x) [/imath] op die massa m uitoefent. De positieve constanten m, [imath] \mathrm{F}_0 [/imath] en [imath] \mathrm{X}_0 [/imath] worden hierbij zo gekozen dat ze voldoen aan: [imath] \mathrm{F}_0 = (\omega_0)^2 \cdot \mathrm{m} \mathrm{X}_0 [/imath] . Op grond van onze eerdere resultaten kiezen we daarbij [imath] \omega_0 = 1.947 \cdot 10^9 \mathrm{rad}/s [/imath]. Dan geldt volgens Newton voor de bewegingen van die massa m dat:

[imath] \mathrm{F}(x) = \mathrm{m} \cdot \ddot{X}(t) [/imath]

[imath] - \mathrm{F}_0 \cdot \arctan(x) = \mathrm{m} \cdot \mathrm{X}_0 \ddot{x}(t) [/imath]

[imath] - \mathrm{F}_0 \cdot \arctan(x) = \mathrm{m} \mathrm{X}_0 \cdot \ddot{x}(t) [/imath]

[imath] - (\omega_0)^2 \cdot \arctan(x) = \ddot{x}(t) [/imath]

[imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x) = 0 [/imath]

Hiermee hebben we dus een mechanisch systeem dat precies voldoet aan onze differentiaalvergelijking [imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x) = 0 [/imath] . Nu is het zaak te bewijzen dat de massa m in dit systeem enkel periodieke bewegingen met een vaste frequentie kan uitvoeren. In dat geval geldt immers vanzelf ook dat de differentiaalvergelijking [imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x) = 0 [/imath] enkel periodieke oplossingen heeft.

 

[ProgHead]

De kracht [imath] \mathrm{F}(x) [/imath] is steeds tegengesteld gericht aan de uitwijking X(t), dus de massa m kan behalve op de positie X=0 nooit blijvend in rust zijn. (Als de massa m op X=0 blijvend in rust is dan is de "beweging" in triviale zin periodiek.) Beweegt de massa m langs de positieve x-as in positieve richting dan ondervindt m een tegenwerkende kracht die de kinetische energie van m gaandeweg uitput. Ver weg van het punt X=0 zal die tegenwerkende kracht bij benadering [imath] - \frac{\pi}{2} \mathrm{F}_0 [/imath] zijn. Voor alle denkbare kinetische energieën van m komt er dus voor m bewegend in positieve richting langs de positieve x-as een punt waarbij m rechtsomkeert moet maken en weer terug in de richting van X=0 moet gaan bewegen. Daarbij zal m dan weer kinetische energie opbouwen omdat de kracht [imath] \mathrm{F}(x) [/imath] bij terug beweging van m in de richting van X=0 die beweging ondersteunt. Aangekomen bij X=0 zal m daar dan ook voorbij schieten. Voor het tegengestelde geval dat de massa m langs de negatieve x-as in negatieve richting beweegt geldt op analoge wijze dat er eveneens in alle gevallen een punt komt waarbij m rechtsomkeert moet maken en weer terug in de richting van X=0 moet gaan bewegen. En ook dan zal m bij X=0 aangekomen daar vanwege de inmiddels weer opgebouwde kinetische energie voorbij schieten.

Kwalitatief is het dus wel duidelijk dat de massa m (behalve in het geval waarin m op X=0 in rust verkeert) om X=0 heen en weer slingert. Nu moeten we alleen nog bewijzen dat bij gegeven beginvoorwaarden iedere daarna volgende slingerbeweging precies hetzelfde verloopt.

 

[ProgHead]

De arbeid die de kracht [imath] \mathrm{F}(x) [/imath] bij beweging van [imath] X = \mathrm{X}_1 [/imath] naar [imath] X = \mathrm{X}_2 [/imath] op de massa m verricht wordt aan de kinetische energie [imath] \mathrm{E}_k [/imath] van die massa m toegevoegd. Dus:

[imath] \Delta \mathrm{E}_k = \int_{ \mathrm{X}_1 }^{ \mathrm{X}_2 } \mathrm{F}(x)\, \mathrm{d} X [/imath]

Het maakt daarbij niet uit of de massa m rechtstreeks van [imath] X = \mathrm{X}_1 [/imath] naar [imath] X = \mathrm{X}_2 [/imath] beweegt of dat m dat via een buiten [imath] \mathrm{X}_1 [/imath] en [imath] \mathrm{X}_2 [/imath] gelegen punt [imath] X = \mathrm{X}_3 [/imath] doet. Immers de kracht [imath] \mathrm{F}(x) [/imath] is voor alle waarden van x gedefinieerd en bovendien naar X en dus ook naar x integreerbaar, zodat:

[imath] \int_{ \mathrm{X}_1 }^{ \mathrm{X}_3 } \mathrm{F}(x)\, \mathrm{d} X \,\,\, + \,\,\, \int_{ \mathrm{X}_3 }^{ \mathrm{X}_2 } \mathrm{F}(x)\, \mathrm{d} X \,\, = \,\, \int_{ \mathrm{X}_1 }^{ \mathrm{X}_2 } \mathrm{F}(x)\, \mathrm{d} X [/imath]

Bovendien hebben we:

[imath] \int_{ \mathrm{X}_1 }^{ \mathrm{X}_1 } \mathrm{F}(x)\, \mathrm{d} X \,\, = \,\, 0 [/imath]

De totale verandering van kinetische energie die de massa m ondergaat wanneer zij na een of meerdere omwegen op een punt [imath] X = \mathrm{X}_1 [/imath] terugkeert is dus nul. En dat geldt voor alle denkbare punten [imath] X = \mathrm{X}_1 [/imath]. Anders gezegd: de kinetische energie van de massa m is van de beginsituatie af enkel nog maar afhankelijk van de positie van de massa m op de x-as en dat specifieke verband verandert niet met de tijd die er vanaf de beginsituatie is verlopen.

We zijn er nu bijna, en in de volgende aflevering hoop ik het bewijs te voltooien.

 

[ProgHead]

We stelden al vast dat de kinetische energie van de massa m van de beginsituatie af enkel nog maar afhankelijk is van de positie van die massa m op de x-as en dat dat specifieke verband ook niet meer verandert met de tijd die er vanaf de beginsituatie is verlopen. Verder was het ook al duidelijk dat de massa m (behalve in het geval waarin m op X=0 in rust verkeert) om X=0 heen en weer slingert. In het geval waarin m op het punt X=0 in rust verkeert volvoert m in triviale zin een periodieke beweging. We hoeven nu dus enkel nog die gevallen te bekijken waarbij m niet op het punt X=0 in rust verkeert. Voor die laatste gevallen ligt de grootte van de snelheid van m voor ieder punt van de x-as dat m aandoet vast via de kinetische energie die bij dat punt hoort. Als beginsituatie dienen zowel de beginsnelheid als de beginpositie van m gegeven te zijn. De bewegingsrichting aan het begin is daarmee bekend, en dus weten we dat m ook in die richting zal blijven bewegen tot het punt waarop de kinetisch energie van m nul wordt. Op dat laatst genoemde punt keert m van bewegingsrichting om en dan zal m in die nieuwe richting blijven bewegen (daarbij eerst weer kinetische energie opbouwende en dan voorbij X=0 weer verliezende) tot de kinetisch energie van m opnieuw nul is. Dan keert m weer van bewegingsrichting om, etc. Een bijzonder geval vormen zulke beginsituaties waarbij de beginpositie weliswaar niet op X=0 ligt maar de beginsnelheid wel nul is. In dat geval start de massa m met een kinetische energie nul en dus op een van de twee mogelijke extreme posities die m voor zo'n beginsituatie op de x-as kan aandoen. Maar voor de rest geldt ook dan hetzelfde verhaal als voor de andere situaties. Zowel de grootte als de richting van de bewegingssnelheid van m liggen dus voor ieder punt dat m heen en weer slingerend op de x-as op zeker moment aandoet vast zodra de beginsituatie gegeven is. Dit laatste kan enkel maar zo zijn als de beweging van m van de beginsituatie af volstrekt periodiek is.

 

[ProgHead]

Hoog tijd om nu te bekijken hoe we het digitale riedeltjes artefact als gewenste feature kunnen gebruiken en aansturen. Digitale artefacten ontstaan (programmeerfouten, onvoldoende RAM, etc. daargelaten) in principe door een in verhouding tot de vereiste precisie onvoldoende bitdiepte en/of te grote sample tijd. Door die twee dingen tijdens het runnen van je programma te regelen zou je de creatie van digitale artefacten moeten kunnen sturen. De vraag is of de online Faust IDE het toelaat om daaraan te sleutelen. Zo niet - dan moet dat met kunst en vliegwerk in het Faust programmaatje zelf gebeuren.

 

[ProgHead]

Heb het progje als test maar even volgestopt met bit crushers maar dat maakt voor de riedeltjes geen enkel verschil!

declare name "crusherRiedel";
import("maths.lib");
import("filters.lib");
import("basics.lib");
start  = 0.01*button("start");
nbits = hslider("resolution",8,1,16,1);
A = _ , start :> bitcrusher(nbits): _*-1 :bitcrusher(nbits): _ ;
B = _ :bitcrusher(nbits): atan : integrator :bitcrusher(nbits): integrator :bitcrusher(nbits): _*1645333507 : bitcrusher(nbits) : _ ;
C = A~B ;
safe(x) = atan((10^-11)*x)/3;
differentiator = zero(1)*SR/1000;
process = C <: _,differentiator : safe,safe : _,_ ;

 

[Disharmonic]

Net geprobeerd, dat is toch raar. Je zou verwachten dat het echt iets anders moet worden met de bitcrusher.

 

[ProgHead]

Ja en ik houd nu mijn hart vast voor wat er (niet) gebeurt bij een andere sample rate..

 

[ProgHead]

Dit geeft je een zekere mate van controle:

declare name "ControlRiedel";

import("maths.lib");
import("filters.lib");

start  = 0.01*button("start");
contr = hslider("control",1,0.0001,30,0.01);

A = _ , start :> _*-1 : _ ;
B = _ : _*contr : atan : _*contr : integrator : _*contr : integrator : _*1645333507 : _ ;
C = A~B ;

safe(x) = atan((10^-11)*x)/3;
differentiator = zero(1)*SR/1000;

process = C <: _,differentiator : safe,safe : _,_ ;

 

[ProgHead]

Zo ziet het fase diagram eruit:

Of zo'n bibberend rechthoekje:

 

[ProgHead]

Is onze differentiaalvergelijking [imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x) = 0 [/imath] wellicht een stijve differentiaalvergelijking?

Stiff equation - Wikipedia

en.wikipedia.org

Dan zouden we mogelijk zulke grappen als de verkregen riedeltjes kunnen verwachten...

 

[Disharmonic]

In het algemeen is numerieke instabiliteit misschien iets om naar te kijken voor onverwachte resultaten.

 

[ProgHead]

Chaos, stijve vergelijkingen, numerieke instabiliteit. Dit zijn hier denk ik verschillende benamingen voor in wezen hetzelfde verschijnsel. Dat de differentiaalvergelijking [imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x) = 0 [/imath] keurige periodieke oplossingen heeft bewijst nog niet dat de daarmee corresponderende differentievergelijking (als discrete pendant van de differentiaalvergelijking) ook keurige periodieke oplossingen heeft. Computers werkten bij het gebruik van numerieke methoden strikt genomen ook niet met differentiaalvergelijkingen maar met de daarmee corresponderende differentievergelijkingen, en het ziet er nu naar uit dat de met [imath] \ddot{x}(t) + (\omega_0)^2 \cdot \arctan(x) = 0 [/imath] corresponderende differentievergelijking wel chaotisch oplossingen heeft. In dat geval doen mijn computer en mijn Faust programma ook niets verkeerd.

 

[ProgHead]

Als we de differentievergelijking weten waarop Faust de simulatie van onze differentiaalvergelijking baseert dan kunnen we de computer ook een stukje van een oplossing van die differentievergelijking laten uitrekenen. Zitten in die oplossing dan ook riedeltjes dan weten we dat die riedeltjes uit het gebruik van een differentievergelijking (in plaats van de differentiaalvergelijking) voortkomen.

 

[ProgHead]

Even terug naar de basis: zie hier de blokken van ons eenvoudige Faust riedeltjes progje:

Houd er ook rekening mee dat de feedback loop een impliciete one-sample delay bevat!

 

[ProgHead]

De "(fi.)integrator" is "pole(1)" dus voldoet aan y(n) = x(n) + y(n-1). Voor de plaats van signalen a, b, c, d en e in B zie de afbeelding hieronder:

We lezen daaruit af dat:

b(n) = b(n-1) + a(n)
c(n) = b(n)/SR
d(n) = d(n-1) + c(n)
e(n) = d(n)/SR

Zodat:

e(n) = (d(n-1)+ c(n))/SR
e(n) = (e(n-1)*SR + c(n))/SR
e(n) = (e(n-1)*SR + b(n)/SR)/SR
e(n) = e(n-1)+ b(n)/SR^2
e(n) = e(n-1)+ (b(n-1) + a(n))/SR^2
e(n) = e(n-1)+ (c(n-1)*SR + a(n))/SR^2
e(n) = e(n-1)+ c(n-1)/SR + a(n)/SR^2
e(n) = e(n-1)+ {d(n-1) - d(n-2)}/SR + a(n)/SR^2
e(n) = e(n-1)+ { e(n-1)*SR - e(n-2)*SR }/SR + a(n)/SR^2
e(n) = e(n-1)+ e(n-1) - e(n-2) + a(n)/SR^2
e(n) = 2e(n-1) - e(n-2) + a(n)/SR^2