solveODE((t,y)=>y, [0,4], 1, {method:'rk23'})
viewof method = Inputs.select(["RK1", "RK2", "RK3", "RK4", "RK12", "RK23", "RK45", "BDF2"], { value:"RK45", label: "Method"})
viewof maxStep = Inputs.range([0.01, 1], { value: 0.5, label: "Step size (or maximum step for variable step methods)", step: 0.01 })
viewof tol = Inputs.range([1e-4, 1], { value: 1e-3, label: "Tolerance", step: 1e-4 })
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})
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 = solveODE(fvdp, [0, 20], [y_start, y_prime_start], { tol, maxStep, method, step:maxStep })
solveODE((t,y)=>y, [0, 4], 1, {method:"rk45", firstStep:-1, tol:1e-3})
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 bdf2(f, T, y0, opt) {
// BDF2 implemented by chatgpt (it doesn't work yet)
const t0 = T[0];
const tf = T[1];
const tolerance = opt.tol ? opt.tol : 1e-3
const h = opt.step ? opt.step : 0.1
const maxIterations = opt.maxIter ? opt.maxIter : 1_000; // Maximum number of iterations
let t = [t0];
let y = [y0];
let hActual = h;
for (let i = 0; i < maxIterations; i++) {
if (t >= tf) {
return y;
}
const yPrev = y;
// BDF2 formula
const yNext = math.add(
math.multiply(
math.divide(
math.subtract(
math.multiply(2, y[i]),
y[i]
),
hActual
),
hActual
),
math.multiply(
math.divide(
math.multiply(3, f(t, y[i])),
2
),
math.pow(hActual, 2)
)
);
// Estimate the error
const error = math.max(math.abs(math.subtract(yNext, y)));
// Adjust the step size based on the error
const scaleFactor = math.divide(tolerance, error);
hActual = math.multiply(h, math.pow(scaleFactor, 1 / 3));
// Update the time and state
t += h;
y = yNext;
}
throw new Error('BDF2 failed to converge.');
}
function _rk(butcherTableau) {
// generates an adaptive runge kutta method from it's butcher tableau
const subtract = math.subtract;
const isUnit = math.isUnit;
const isBigNumber = math.isBigNumber;
const isNumeric = math.isNumeric
const bignumber = math.bignumber
const divide = math.divide
const add = math.add
const multiply = math.multiply
const max = math.max
const abs = math.abs
const map = math.map
const smaller = math.smaller
const larger = math.larger
const isPositive = math.isPositive
return function (f, tspan, y0, options) {
// adaptive runge kutta methods
const hasBigNumbers = [
...tspan,
...y0,
options.maxStep,
options.minStep
].some(isBigNumber);
const wrongTSpan = !(
tspan.length === 2 &&
((isNumeric(tspan[0]) && isNumeric(tspan[1])) ||
(isUnit(tspan[0]) && isUnit(tspan[1])))
);
if (wrongTSpan) {
throw new Error(
'"tspan" must be an Array of two numeric values or two units [tStart, tEnd]'
);
}
const t0 = tspan[0]; // initial time
const tf = tspan[1]; // final time
const steps = 1; // divide 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
const [a, c, b, bp] = hasBigNumbers
? [
bignumber(butcherTableau.a),
bignumber(butcherTableau.c),
bignumber(butcherTableau.b),
bignumber(butcherTableau.bp)
]
: [
butcherTableau.a,
butcherTableau.c,
butcherTableau.b,
butcherTableau.bp
];
let h = options.firstStep
? options.firstStep
: divide(subtract(tf, t0), steps); // define the first step size
const t = [t0]; // start the time array
const y = [y0]; // start the solution array
const dletaB = 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)) {
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(add(t[n], multiply(c[i], h)), add(y[n], multiply(h, a[i], k)))
);
}
// estimate the error by comparing solutions of different orders
const TE = max(
abs(map(multiply(dletaB, k), (X) => (isUnit(X) ? X.value : X)))
);
if (TE < tol && tol / TE > 1 / 4) {
// if within tol then push everything
t.push(add(t[n], h));
y.push(add(y[n], multiply(h, b, k)));
n++;
}
// estimate the delta value that will affect the step size
let delta = 0.84 * (tol / TE) ** (1 / 5);
if (smaller(delta, minDelta)) {
delta = minDelta;
} else if (larger(delta, maxDelta)) {
delta = maxDelta;
}
delta = hasBigNumbers ? bignumber(delta) : delta;
h = multiply(h, delta);
if (maxStep && larger(abs(h), abs(maxStep))) {
h = isPositive(h) ? abs(maxStep) : multiply(-1, abs(maxStep));
} else if (minStep && smaller(abs(h), abs(minStep))) {
h = isPositive(h) ? abs(minStep) : multiply(-1, abs(minStep));
}
iter++;
if (iter > maxIter) {
throw new Error(
"Maximum number of iterations reached, try changing options"
);
}
}
return { t, y };
};
}
function _rkdp(f, T, y0, options) {
// Dormand Prince method
// Define the butcher table
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];
const butcherTableau = {a, c, b, bp}
// Solve an adaptive step size rk method
return _rk(butcherTableau)(f, T, y0, options)
}
function rk23(f, T, y0, options) {
// Bogacki–Shampine method
// Define the butcher table
const a = [
[],
[1 / 2],
[0, 3 / 4],
[2 / 9, 1 / 3, 4 / 9]
];
const c = [null, 1 / 2, 3 / 4, 1];
const b = [2/9, 1/3, 4/9, 0];
const bp = [7/24, 1/4, 1/3, 1/8];
const butcherTableau = {a, c, b, bp}
// Solve an adaptive step size rk method
return _rk(butcherTableau)(f, T, y0, options)
}
function rk12(f, T, y0, options) {
// Fehlberg rk12
// Define the butcher table
const a = [
[],
[1 / 2],
[1 / 256, 255 / 256]
];
const c = [null, 1 / 2, 1];
const b = [1/512, 255/256, 1/512];
const bp = [1/256, 255/256, 0];
const butcherTableau = {a, c, b, bp}
// Solve an adaptive step size rk method
return _rk(butcherTableau)(f, T, y0, options)
}
function euler(f, T, y0, opt) {
// https://mathworld.wolfram.com/EulerForwardMethod.html
const steps = opt.steps ? opt.steps : 10;
const t0 = T[0]; // initial time
const tf = T[1]; // final time
let h = opt.step ? opt.step : math.divide(math.subtract(tf, t0), steps);
let y = [y0];
let t = [t0];
for (let n = 0; ongoing(t[n], tf, h); n++) {
// trim the time step so that it doesn't overshoot
h = trimStep(t[n], tf, h);
t.push(math.add(t[n], h));
y.push(math.add(y[n], math.multiply(h, f(t[n], y[n]))));
}
return { t, y };
}
solveODE((t,y)=>y, [0,4], [1, 1-1e-6], {method:'RK3', step:1})
function midpoint(f, T, y0, opt) {
// https://mathworld.wolfram.com/Runge-KuttaMethod.html
const steps = opt.steps ? opt.steps : 10;
const t0 = T[0]; // initial time
const tf = T[1]; // final time
let h = opt.step ? opt.step : math.divide(math.subtract(tf, t0), steps);
let y = [y0];
let t = [t0];
for (let n = 0; ongoing(t[n], tf, h); n++) {
// trim the time step so that it doesn't overshoot
let k = [];
h = trimStep(t[n], tf, h);
k.push(math.dotMultiply(h, f(t[n], y[n])))
k.push(math.dotMultiply(h, f(math.add(t[n], math.divide(h, 2)), math.add(
y[n],
math.divide(
k[0],
2
)))))
t.push(math.add(t[n], h));
y.push(math.add(y[n], k[1]));
}
return { t, y };
}
function RK4(f, T, y0, opt) {
// https://mathworld.wolfram.com/Runge-KuttaMethod.html
const steps = opt.steps ? opt.steps : 10;
const t0 = T[0]; // initial time
const tf = T[1]; // final time
const b = [1/6, 1/3, 1/3, 1/6]
let h = opt.step ? opt.step : math.divide(math.subtract(tf, t0), steps);
let y = [y0];
let t = [t0];
for (let n = 0; ongoing(t[n], tf, h); n++) {
let k = []
h = trimStep(t[n], tf, h);
k.push(math.multiply(h, f(t[n], y[n])))
k.push(math.multiply(h, f(
math.add(t[n], math.divide(h, 2)),
math.add(y[n], math.divide(k[0], 2))
)));
k.push(math.multiply(h, f(math.add(t[n], math.divide(h, 2)), math.add(
y[n],
math.divide(k[1], 2)
))));
k.push(math.multiply(h, f(math.add(t[n], h), math.add(y[n], k[2]))));
t.push(math.add(t[n], h));
y.push(math.add(y[n], math.multiply(b, k)));
}
return { t, y }
}
function RK3(f, T, y0, opt) {
// https://mathworld.wolfram.com/Runge-KuttaMethod.html
const steps = opt.steps ? opt.steps : 10;
const t0 = T[0]; // initial time
const tf = T[1]; // final time
const b = [2/9, 1/3, 4/9]
let h = opt.step ? opt.step : math.divide(math.subtract(tf, t0), steps);
let y = [y0];
let t = [t0];
for (let n = 0; ongoing(t[n], tf, h); n++) {
let k = []
h = trimStep(t[n], tf, h);
k.push(math.multiply(h, f(t[n], y[n])))
k.push(math.multiply(h, f(
math.add(t[n], math.divide(h, 2)),
math.add(y[n], math.divide(k[0], 2))
)));
k.push(math.multiply(h, f(math.add(t[n], math.multiply(h, 3/4)), math.add(
y[n],
math.multiply(k[1], 3/4)
))));
t.push(math.add(t[n], h));
y.push(math.add(y[n], math.multiply(b, k)));
}
return { t, y }
}
function _solveODE(f, T, y0, opt) {
const method = opt.method ? opt.method : "RK45";
const methods = {
RK1: euler,
RK2: midpoint,
RK3,
RK4,
RK45: _rkdp,
AB1: euler,
RK23: rk23,
RK12: rk12,
BDF2: bdf2
};
let methodOptions = { ...opt }; // clone the options object
delete methodOptions.method; // delete the method as it won't be needed
return methods[method.toUpperCase()](f, T, y0, methodOptions);
}
solveODE = math.typed("solveODE", {
"function, Array, Array, Object": _solveODE,
"function, Array, Array": (f, T, y0) => _solveODE(f, T, y0, {}),
"function, Array, number | BigNumber": (f, T, y0) => {
const sol = _solveODE(f, T, [y0], {});
return { t: sol.t, y: sol.y.map((Y) => Y[0]) };
},
"function, Array, number | BigNumber, Object": (f, T, y0, options) => {
const sol = _solveODE(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(t_n, tf, h) {
return math.isPositive(h) ? math.smaller(t_n, tf) : math.larger(t_n, tf);
}
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.
