Hello,
I implemented it yesterday "just for fun" (it took me the whole day!)
and surprisingly the code doesn't look bad to me. It is very concise
and easy to extend.
Of course this is just a toy, but perhaps it can be useful to someone
so I decided to share it.
Usage: rk(order, n-dimensional-system, (x0,x1, ..., xn));
Simple example, spring-mass oscillator:
M * x''(t) = -K * x(t);
As usual, let v(t) = x'(t)
x'(t) = v(t);
v'(t) = -K / M * x(t);
the faust code:
rk = library("rk.dsp");
spring(M,K, t,x,v) = v, -K / M * x;
process = rk.rk(
4, // order: 1-4
spring(10,40), // parameters: M = 10, K = 40
(1, 0) // initial values: x(0) = 1, v(0) = 0
);
Note that rk() accepts the differential system of any dimension.
Note also that rk.dsp implements the single generic operator, rk_step()
which takes the Butcher tableau in arguments. So if you want, for example,
the 3/8-rule fourth-order method described by
0 |
1/3 | 1/3
2/3 |-1/3 1
1 | 1 -1 1
----------------
| 1/8 3/8 3/8 1/8
tableau you can just use
rk_step(
// 1st arg: c1 .. c3
(0, 1/3, 2/3, 1),
// 2nd arg: a21 .. a43
( 1/3,
-1/3, 1,
1, -1, 1,
// and b1 .. b4 to avoid another argument
1/8, 3/8, 3/8, 1/8)
);
Oleg.
-------------------------------------------------------------------------------
import("stdfaust.lib");
mapl(op, l) = par(i, ba.count(l), op(ba.take(i+1, l)));
loop(op, ini) = op ~ mapl(i, ini) with {
i(0) = _; // optimization
i(x) = +(prefix(x,0));
};
rk_step(ts,ks, t,h,op) = si.bus(N)
: seq(i, O, pushk(i+1, vec(i), opt(i)))
: ro.interleave(N,O+1) : par(i,N, vec(O))
with {
O = ba.count(ts);
N = inputs(op)-1;
pushk(K, vec, op) = si.bus(K*N) <: si.bus(K*N), vecop
with {
vecop = ro.interleave(N,K) : par(i,N, vec) : op;
};
vec(0) = _;
vec(k) = _, mapl(*, ba.subseq(ks, k*(k-1)/2, k)) :> _;
opt(k) = op(t + h * ba.take(k+1,ts)) : par(i,N, *(h));
};
rk_th(O, t,h,op) = loop(step(O, t,h,op)) with {
// pass Butcher tableau in (ts, ks)
step(1) = rk_step((0), (1));
step(2) = rk_step((0,1/2), (1/2, 0,1));
step(3) = rk_step((0,1/2,1), (1/2,-1,2, 1/6,2/3,1/6));
step(4) = rk_step((0,1/2,1/2,1), (1/2,0,1/2,0,0,1, 1/6,1/3,1/3,1/6));
};
rk(O) = rk_th(O, ba.time/ma.SR, 1/ma.SR);
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Faudiostream-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/faudiostream-users
