Public
Edited
May 2, 2023
Insert cell
Insert cell
function timeRange(t0, tf, h){
const minSteps = 8
let t = [t0]
// if no time step was supplied make one with the time range
h = h ? h : math.divide(math.subtract(tf, t0), minSteps)
let n = 0
while(notFinished(t[n], tf, h)){
h = trimStep(t[n], tf, h)
t.push(math.add(t[n], h))
n ++
}
// return an empty array if it didn't iterate a single time
t = t.length > 1 ? t : []
return t
}
Insert cell
odeEuler = function (f, T, y0) {
// https://mathworld.wolfram.com/EulerForwardMethod.html
const t = timeRange(...T)
const N = t.length - 1 // number of times the method has to run

let y = [y0]

for (let n = 0; n < N; n++) {
const h = math.subtract(t[n + 1], t[n])
y.push(
math.add(
y[n],
math.dotMultiply(
h,
f(t[n], y[n])
)
)
)
}
return { t, y }
}
Insert cell
odeEuler((t,y)=>y, [0,4,1], [1]).y
Insert cell
odeEuler((t,y) => y, [0,4,1], 1).y
Insert cell
function notFinished(t_n, tf, h){
// Checks if the time has reached the final time in the direction of the time step
if(math.isPositive(h)){
// If time step is positive, return true if
return math.smaller(t_n, tf)
}
else if(math.isNegative(h)){
return math.smaller(tf, t_n)
}
}
Insert cell
function trimStep(t_n, tf, h){
// Trims the next step to reach tf (not after tf)
const next = math.add(t_n,h)
if(math.isPositive(h)){
if(math.larger(next, tf)){
h = math.subtract(tf, t_n)
}
}
else if(math.isNegative(h)){
if(math.smaller(next, tf)){
h = math.subtract(tf, t_n)
}
}
return h
}
Insert cell
function ode45(f, T, y0){
// Dormand Prince
const tol = 1e-4
const hmax = 1
const t0 = T[0]
const tf = T[1]
const max_iter = 10_000
//const unit = math.isUnit(y0) ? y0.formatUnits() : null
const minSteps = 8
//const units = math.map(y0, (y)=> math.isUnit(y) ? y.formatUnits() : null)
let h = T[2] ? T[2] : math.divide(math.subtract(T[1],T[0]), minSteps)
let t = [t0]
let y = [y0]

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 = [0, 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 delta_b = math.subtract(b, bp) // b - bp

let n = 0
let iter = 0
while (notFinished(t[n], tf, h) && (iter < max_iter)){
let k = []

h = trimStep(t[n], tf, h)

k.push(f(t[n], y[n]))

k.push(f(math.add(t[n], math.dotMultiply(c[1], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.dotMultiply(a[1][0], k[0])
)
)
))

k.push(f(
math.add(t[n], math.dotMultiply(c[2], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[2][0], k[0]),
math.dotMultiply(a[2][1], k[1])
)
)
)
))

k.push(f(
math.add(t[n], math.dotMultiply(c[3], h))
,
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[3][0], k[0]),
math.dotMultiply(a[3][1], k[1]),
math.dotMultiply(a[3][2], k[2])
)
)
)
))


k.push(f(
math.add(t[n], math.dotMultiply(c[4], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[4][0], k[0]),
math.dotMultiply(a[4][1], k[1]),
math.dotMultiply(a[4][2], k[2]),
math.dotMultiply(a[4][3], k[3])
)
)
)
))

k.push(f(
math.add(t[n], h),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[5][0], k[0]),
math.dotMultiply(a[5][1], k[1]),
math.dotMultiply(a[5][2], k[2]),
math.dotMultiply(a[5][3], k[3]),
math.dotMultiply(a[5][4], k[4])
)
)
)
))

k.push(f(
math.add(t[n], h),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[6][0], k[0]),
math.dotMultiply(a[6][1], k[1]),
math.dotMultiply(a[6][2], k[2]),
math.dotMultiply(a[6][3], k[3]),
math.dotMultiply(a[6][4], k[4]),
math.dotMultiply(a[6][5], k[5])
)
)
)
))

const TE = math.abs(
math.add(
math.dotMultiply(delta_b[0], k[0]),
math.dotMultiply(delta_b[2], k[2]),
math.dotMultiply(delta_b[3], k[3]),
math.dotMultiply(delta_b[4], k[4]),
math.dotMultiply(delta_b[5], k[5]),
math.dotMultiply(delta_b[6], k[6])
)
)


const TEnum = math.isUnit(TE) ? TE.toNumber(TE.formatUnits()) : TE

const delta =
0.84 *
math.dotPow(
math.dotDivide(tol, TEnum),
(1/5)
)

if ( TEnum < tol ) {
t.push(math.add(t[n], h))
y.push(
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(b, k)
)
)
)
n += 1
}
if (delta <= 0.1){
h = math.multiply(h, 0.1)
}
else if( delta >= 4)
{
h = math.multiply(h, 4)
}
else{
h = math.multiply(h, delta)
}
if ( h > hmax){
h = hmax
}
//h = math.multiply(h, delta)
iter ++
}
return {t,y}
}
Insert cell
({method:'rk45', tol:1e-4})
Insert cell
function ode(f, T, y0){
// Dormand Prince
const tol = 1e-4
const hmax = 1
const t0 = T[0]
const tf = T[1]
const max_iter = 10_000
const minSteps = 8
const units = math.map(y0, (y)=> math.isUnit(y) ? y.formatUnits() : "")
const timeUnits = math.isUnit(t0) ? t0.formatUnits() : ""
const kUnits = math.dotDivide(math.unit(units), math.unit(timeUnits))

let h = T[2] ? T[2] : math.divide(math.subtract(T[1],T[0]), minSteps)
let t = [t0]
let y = [y0]

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 = [0, 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 delta_b = math.subtract(b, bp) // b - bp

let n = 0
let iter = 0
while (notFinished(t[n], tf, h) && (iter < max_iter)){
let k = []

h = trimStep(t[n], tf, h)

k.push(f(t[n], y[n]))

k.push(f(math.add(t[n], math.dotMultiply(c[1], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.dotMultiply(a[1], k[0])
)
)
))

k.push(f(
math.add(t[n], math.dotMultiply(c[2], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(a[2], k)
)
)
))

k.push(f(
math.add(t[n], math.dotMultiply(c[3], h))
,
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(a[3], k)
)
)
))


k.push(f(
math.add(t[n], math.dotMultiply(c[4], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(a[4], k)
)
)
))

k.push(f(
math.add(t[n], h),
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(a[5], k)
)
)
))

k.push(f(
math.add(t[n], h),
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(a[6], k)
)
)
))

const TE = math.abs(
math.multiply(delta_b, k)
)

const TEnum = math.hypot(
TE.map((y, index) => math.isUnit(y) ? y.toNumber(kUnits[index]) : y)
)

const delta =
0.84 *
math.dotPow(
math.dotDivide(tol, TEnum),
(1/5)
)

if ( TEnum < tol ) {
t.push(math.add(t[n], h))
y.push(
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(b, k)
)
)
)
n += 1
}
/*
if (delta <= 0.1){
h = math.multiply(h, 0.1)
}
else if( delta >= 4)
{
h = math.multiply(h, 4)
}
else{
h = math.multiply(h, delta)
}
if ( h > hmax){
h = hmax
}
*/
h = math.multiply(h, delta)
iter ++
}
return {t,y}
}
Insert cell
solArray = ode((t,y) => math.multiply(y, math.unit('1/s')), [math.unit(0,'s'), math.unit(4,'s')], [math.unit("1m")])
Insert cell
solArray1 = ode((t,y) => y, [0, 4], [1,1.001])
Insert cell
math.format(solArray.y[solArray.y.length-1])
Insert cell
sol = ode45((t,y) => y, [0, 4], 1)
Insert cell
solUnit = ode45((t,y) => math.multiply(y, math.unit('1/s')), [math.unit(0,'s'), math.unit(4,'s')], math.unit("1m"))
Insert cell
math.exp(4)
Insert cell
sol.y[sol.y.length-1]
Insert cell
math.format(solUnit.y[sol.y.length-1])
Insert cell
math.exp(4)
Insert cell
math = require('mathjs@11.8')
Insert cell

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.
Learn more