78 lines
2.9 KiB
Lua
78 lines
2.9 KiB
Lua
-- dop853.lua: Dormand-Prince 8th order (13 NFE, fixed step)
|
|
|
|
solver = {
|
|
name = "dop853",
|
|
display = "DOP853 (13 NFE)",
|
|
description = "8th order Dormand-Prince (maximum accuracy)",
|
|
nfe = 13,
|
|
order = 8,
|
|
needs_model = true,
|
|
stateful = false,
|
|
stochastic = false,
|
|
}
|
|
|
|
local C = {
|
|
1/18, 1/12, 1/8, 5/16, 3/8,
|
|
59/400, 93/200, 5490023248/9719169821,
|
|
13/20, 1201146811/1299019798,
|
|
1, 1,
|
|
}
|
|
local A = {
|
|
{1/18},
|
|
{1/48, 1/16},
|
|
{1/32, 0, 3/32},
|
|
{5/16, 0, -75/64, 75/64},
|
|
{3/80, 0, 0, 3/16, 3/20},
|
|
{29443841/614563906, 0, 0, 77736538/692538347, -28693883/1125000000, 23124283/1800000000},
|
|
{16016141/946692911, 0, 0, 61564180/158732637, 22789713/633445777, 545815736/2771057229, -180193667/1043307555},
|
|
{39632708/573591083, 0, 0, -433636366/683701615, -421739975/2616292301, 100302831/723423059, 790204164/839813087, 800635310/3783071287},
|
|
{246121993/1340847787, 0, 0, -37695042795/15268766246, -309121744/1061227803, -12992083/490766935, 6005943493/2108947869, 393006217/1396673457, 123872331/1001029789},
|
|
{-1028468189/846180014, 0, 0, 8478235783/508512852, 1311729495/1432422823, -10304129995/1701304382, -48777925059/3047939560, 15336726248/1032824649, -45442868181/3398467696, 3065993473/597172653},
|
|
{185892177/718116043, 0, 0, -3185094517/667107341, -477755414/1098053517, -703635378/230739211, 5731566787/1027545527, 5232866602/850066563, -4093664535/808688257, 3962137247/1805957418, 65686358/487910083},
|
|
{403863854/491063109, 0, 0, -5068492393/434740067, -411421997/543043805, 652783627/914296604, 11173962825/925320556, -13158990841/6184727034, 3936647629/1978049680, -160528059/685178525, 248638103/1413531060, 0},
|
|
}
|
|
local B = {
|
|
14005451/335480064, 0, 0, 0, 0,
|
|
-59238493/1068277825, 181606767/758867731, 561292985/797845732,
|
|
-1041891430/1371343529, 760417239/1151165299, 118820643/751138087,
|
|
-528747749/2220607170, 1/4,
|
|
}
|
|
|
|
function step(xt, vt, t_curr, t_prev, n, model_fn, vt_buf)
|
|
local dt = t_curr - t_prev
|
|
local xt_orig = {}
|
|
local k1 = {}
|
|
for i = 0, n-1 do xt_orig[i] = xt[i]; k1[i] = vt[i] end
|
|
|
|
local ks = {k1}
|
|
local x_tmp = {}
|
|
|
|
-- 12 extra stages
|
|
for s = 1, 12 do
|
|
local a_row = A[s]
|
|
for i = 0, n-1 do
|
|
local combo = 0
|
|
for j = 1, #a_row do
|
|
if a_row[j] ~= 0 and ks[j] then
|
|
combo = combo + a_row[j] * ks[j][i]
|
|
end
|
|
end
|
|
x_tmp[i] = xt_orig[i] - dt * combo
|
|
end
|
|
-- Write x_tmp to xt for model_fn
|
|
for i = 0, n-1 do xt[i] = x_tmp[i] end
|
|
model_fn(xt, t_curr - C[s] * dt)
|
|
ks[s+1] = {}
|
|
for i = 0, n-1 do ks[s+1][i] = vt_buf[i] end
|
|
end
|
|
|
|
-- Combine
|
|
for i = 0, n-1 do
|
|
local sol = 0
|
|
for j = 1, 13 do
|
|
if B[j] ~= 0 and ks[j] then sol = sol + B[j] * ks[j][i] end
|
|
end
|
|
xt[i] = xt_orig[i] - dt * sol
|
|
end
|
|
end
|