Public
Edited
Jul 5, 2023
1 fork
viewof y_start = Inputs.range([0.1, 2.5], {value:1, label: "Initial y value", step: 0.01})
viewof y_prime_start = Inputs.range([-2.5, 2.5], {value:0, label: "Initial dy/dx value", step: 0.01})
viewof miu = Inputs.range([0.1, 4], {value:1, label: "Miu Value", step: 0.01})
viewof tol = Inputs.range([1e-5, 1e-1], {value:1e-3, label: "Tolerance", step: 1e-4})
fvdp = (t, y) => {
// this can also be defined in one line as the reference
// [y[1], miu*(1 - y[0]**2) * y[1] - y[0] ]
const [x, dxdt] = y
return [dxdt, miu*(1 - x**2) * dxdt - x]
}
sol = rkdp(fvdp, [0, 20], [y_start, y_prime_start], {tol:tol})
vdp = sol.y.map((X, i) => ({t:sol.t[i], x:X[0], dxdt:X[1]}))
// an alternative using arquero
// vdp2 = aq.table({t: sol.t, x:sol.y.map(X=>X[0]), dxdt:sol.y.map(X=>X[1])})
viewof table = Inputs.table(vdp)
function _rkdp(f, T, y0, options) {
// Dormand Prince method
// Define the butcher tableau
const a = [
[],
[1 / 5],
[3 / 40, 9 / 40],
[44 / 45, -56 / 15, 32 / 9],
[19372 / 6561, -25360 / 2187, 64448 / 6561, -212 / 729],
[9017 / 3168, -355 / 33, 46732 / 5247, 49 / 176, -5103 / 18656],
[35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84]
];
const c = [null, 1 / 5, 3 / 10, 4 / 5, 8 / 9, 1, 1];
const b = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84, 0];
const bp = [5179 / 57600, 0, 7571 / 16695, 393 / 640, -92097 / 339200, 187 / 2100, 1 / 40];
// Solve an adaptive step size rk method
return _rk(a, c, b, bp)(f, T, y0, options)
}
math.typeOf(math.unit('1m'))
function _rk(a, c, b, bp) {
// generates an adaptive runge kutta method from it's butcher tableau
return function (f, T, y0, options){
// adaptive runge kutta methods
const t0 = T[0]; // initial time
const tf = T[1]; // final time
const steps = 1; // divide the time in this number of steps
const tol = options.tol ? options.tol : 1e-4; // define a tolerance (must be an option)
const maxStep = options.maxStep;
const minStep = options.minStep;
const minDelta = options.minDelta ? options.minDelta : 0.2;
const maxDelta = options.maxDelta ? options.maxDelta : 5;
const maxIter = options.maxIter ? options.maxIter : 10_000; // stop inifite evaluation if something goes wrong
let h = options.firstStep ?
options.firstStep :
math.divide(math.subtract(tf, t0), steps); // define the initial step size
const t = [t0]; // start the time array
const y = [y0]; // start the solution array
const delta_b = math.subtract(b, bp); // b - bp
let n = 0;
let iter = 0;
// iterate unitil it reaches either the final time or maximum iterations
while (ongoing(t[n], tf, h) && iter < maxIter) {
const k = [];
// trim the time step so that it doesn't overshoot
h = trimStep(t[n], tf, h);
// calculate the first value of k
k.push(f(t[n], y[n]));
// calculate the rest of the values of k
for (let i = 1; i < c.length; ++i) {
k.push(
f(
math.add(t[n], math.multiply(c[i], h)),
math.add(y[n], math.multiply(h, a[i], k))
)
);
}
// estimate the error by comparing solutions of different orders
const TE = math.max(
math.abs(
math.map(math.multiply(delta_b, k), (X) =>
math.isUnit(X) ? X.value : X
)
)
);
if (TE < tol && tol / TE > 1 / 4) {
// if within tol then push everything
t.push(math.add(t[n], h));
y.push(math.add(y[n], math.multiply(h, b, k)));
n ++;
}
// estimate the delta value that will affect the step size
const delta = 0.84 * (tol / TE) ** (1 / 5);
if (delta < minDelta) {
h = math.multiply(h, minDelta);
} else if (delta > maxDelta) {
h = math.multiply(h, maxDelta);
} else {
h = math.multiply(h, delta);
}
if (maxStep && math.larger(math.abs(h), math.abs(maxStep))) {
h = math.isPositive(h) ? math.abs(maxStep) : math.multiply(-1, math.abs(maxStep));
} else if (minStep && math.smaller(math.abs(h), math.abs(minStep))) {
h = math.isPositive(h) ? math.abs(minStep) : math.multiply(-1, math.abs(minStep));
}
iter++;
}
return { t, y };
}
}
rkdp = math.typed('rkdp',{
"function, Array, Array, Object": _rkdp,
"function, Array, Array": (f, T, y0) => _rkdp(f, T, y0, {}),
"function, Array, number | BigNumber | Unit": (f, T, y0) => {
const sol = _rkdp(f, T, [y0], {})
return {t:sol.t, y: sol.y.map(Y => Y[0])}
},
"function, Array, number | BigNumber | Unit, Object": (f, T, y0, options) => {
const sol = _rkdp(f, T, [y0], options)
return {t:sol.t, y: sol.y.map(Y => Y[0])}
}
})
// returns true if the time has not reached tf for both postitive an negative step (h)
function ongoing(tn, tf, h) {
return math.isPositive(h) ?
math.smaller(tn, tf) :
math.larger(tn, tf)
}
rkdp((t,y)=>y, [0,4],[1], {firstStep:1})
function trimStep(t_n, tf, h){
// Trims the time step so that the next step doesn't overshoot
const next = math.add(t_n, h)
return (
(math.isPositive(h) && math.larger(next, tf)) ||
(math.isNegative(h) && math.smaller(next, tf))
) ?
math.subtract(tf, t_n) :
h
}
math = require('mathjs@11.8')
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
