SourceBrowser

meta-imp-functions.lmt /size: 8647 b last modification: 2024-01-16 10:22
1if not modules then modules = { } end modules ['meta-imp-functions'] = {
2    version   = 1.001,
3    comment   = "companion to meta-imp-functions.mkxl",
4    author    = "Hans Hagen, PRAGMA-ADE, Hasselt NL",
5    copyright = "PRAGMA ADE / ConTeXt Development Team",
6    license   = "see context related readme files",
7}
8
9local formatters = string.formatters
10local sequenced  = table.sequenced
11
12local noffunctions = 0
13local version      = 1
14
15local function preparecache(p)
16    noffunctions = noffunctions + 1
17local action = p.action
18p.action = nil
19    local hash = md5.HEX(sequenced(p))
20    local name = formatters["mkiv-%s-m-f-%03i.lua"](tex.jobname,noffunctions)
21p.action = action
22    return name, hash
23end
24
25local function getcache(p)
26    local cache = p.cache
27    if cache then
28        local name, hash = preparecache(p)
29        local data = table.load(name)
30        if data and data.hash == hash and data.version == version and data.data then
31            return hash, name, data.data
32        else
33            return hash, name, false
34        end
35    else
36        return false, false, false
37    end
38end
39
40local function setcache(hash,name,data)
41    local result = {
42        version = version,
43        hash    = hash,
44        data    = data,
45    }
46    table.save(name,result)
47end
48
49local injectpath      = mp.inject.path
50local getparameterset = metapost.getparameterset
51
52local report = logs.reporter("metapost","functions")
53
54local functions = { }
55local actions   = { }
56
57function mp.registerfunction(specification)
58    local name   = specification.name
59    functions[name] = specification
60end
61
62function mp.registeraction(specification)
63    local name   = specification.name
64    actions[name] = specification
65end
66
67metapost.registerscript("processfunction", function()
68    local specification = getparameterset("function")
69    local name    = specification.name
70    local lua     = specification.lua
71    local fnction = functions[name]
72    local action  = lua and actions[lua]
73    if fnction then
74        if action then
75            specification.action = action.action
76        end
77     -- statistics.starttiming(functions)
78        fnction.action(specification)
79     -- statistics.stoptiming(functions)
80    end
81end)
82
83-- statistics.register("mp function time", function()
84--     return statistics.elapsedseconds(functions,"including feedback to metapost")
85-- end)
86
87-- Here comes the fancy stuff:
88
89local math = math
90local sqrt = math.sqrt
91
92local mathfunctions = math.functions or { }
93math.functions      = mathfunctions
94
95-- Todo : reference where we got the factors from because those from
96--
97-- This is Runge-Kutta-Merson 4("5")
98-- See Table 4.1. Merson 4("5") of Hairer, Nørsett, Wanner - Solving Ordinary Differential Equations I (Springer, 2008)
99--
100-- function mathfunctions.rungekutta(specification)
101--     local f      = specification.action or function(t,x,y) return x, y end
102--     local x      = specification.x      or 0
103--     local y      = specification.y      or 0
104--     local t      = 0
105--     local tmax   = specification.tmax   or 1
106--     local dt     = specification.dt     or tmax/10
107--     local eps    = specification.eps    or dt/10
108--     local r      = 1
109--  -- local result = { { x, y, x, y, x, y } }
110--     local result = { { x, y } }
111--     while t < tmax do
112--         local k1x, k1y = f(t,              x,
113--                                            y)
114--         k1x = dt * k1x
115--         k1y = dt * k1y
116--         local k2x, k2y = f(t + (1/3) * dt, x + (1/3) * k1x,
117--                                            y + (1/3) * k1y)
118--         k2x = dt * k2x
119--         k2y = dt * k2y
120--         local k3x, k3y = f(t + (1/3) * dt, x + (1/6) * k1x + (1/6) * k2x,
121--                                            y + (1/6) * k1y + (1/6) * k2y)
122--         k3x = dt * k3x
123--         k3y = dt * k3y
124--         local k4x, k4y = f(t + (1/2) * dt, x + (1/8) * k1x + (3/8) * k3x,
125--                                            y + (1/8) * k1y + (3/8) * k3y)
126--         k4x = dt * k4x
127--         k4y = dt * k4y
128--         local k5x, k5y = f(t +         dt, x + (1/2) * k1x - (3/2) * k3x - (2) * k4x,
129--                                            y + (1/2) * k1y - (3/2) * k3y - (2) * k4y)
130--         k5x = dt * k5x
131--         k5y = dt * k5y
132--         --
133--         local teps = sqrt(((1/10-1/6) * k1x + (3/10) * k3x + (2/5-2/3) * k4x + (1/5 -1/6) * k5x)^2 +
134--                           ((1/10-1/6) * k1y + (3/10) * k3y + (2/5-2/3) * k4y + (1/5 -1/6) * k5y)^2 )
135--         if teps < eps then
136--             dt = 0.9 * dt * (eps/teps)^(1/4)
137--             x = x + (1/10) * k1x + (3/10) * k3x + (2/5) * k4x + (1/5) * k5x
138--             y = y + (1/10) * k1y + (3/10) * k3y + (2/5) * k4y + (1/5) * k5y
139--             r = r + 1
140--          -- result[r] = { x, y, x, y, x, y }
141--             result[r] = { x, y }
142--             t = t + dt
143--         else
144--             dt = 0.9 * dt * (eps/teps)^(1/3)
145--         end
146--     end
147--     return result
148-- end
149
150local function rungekutta(specification)
151    local f        = specification.action   or function(t,x,y) return x, y end
152    local x        = specification.x        or 0
153    local y        = specification.y        or 0
154    local tmin     = specification.tmin     or 0
155    local tmax     = specification.tmax     or 1
156    local t        = tmin
157    local rmax     = specification.maxpath  or 0
158    local stepsize = specification.stepsize or "adaptive"
159    local dt       = specification.dt       or (tmax-tmin)/10
160    local eps      = specification.eps      or dt/10
161    local kind     = specification.kind     or specification.type    -- xy x y
162    local adaptive = stepsize == "adaptive"
163    local r        = 1
164    local result
165    if kind ~= "tx" and kind ~= "ty" then
166        kind = "xy"
167    end
168    if kind == "xy" then
169     -- result = { { x, y, x, y, x, y } }
170        result = { { x, y } }
171    elseif kind == "tx" then
172     -- result = { { x, x, t, x, t, x } }
173        result = { { t, x } }
174    else
175     -- result = { { x, y, t, y, t, y } }
176        result = { { t, y } }
177    end
178    local hash, name, data = getcache(specification)
179    if data then
180     -- print(hash,name,"REUSING")
181        return data
182    else
183     -- print(hash,name,"GENERATING")
184    end
185    if rmax == 0 then
186        rmax = 0xFFFF
187    end
188
189    while t < tmax do
190        local k1x, k1y = f(t,              x,
191                                           y)
192        k1x = dt * k1x
193        k1y = dt * k1y
194        local k2x, k2y = f(t + (1/3) * dt, x + (1/3) * k1x,
195                                           y + (1/3) * k1y)
196        k2x = dt * k2x
197        k2y = dt * k2y
198        local k3x, k3y = f(t + (1/3) * dt, x + (1/6) * k1x + (1/6) * k2x,
199                                           y + (1/6) * k1y + (1/6) * k2y)
200        k3x = dt * k3x
201        k3y = dt * k3y
202        local k4x, k4y = f(t + (1/2) * dt, x + (1/8) * k1x + (3/8) * k3x,
203                                           y + (1/8) * k1y + (3/8) * k3y)
204        k4x = dt * k4x
205        k4y = dt * k4y
206        local k5x, k5y = f(t +         dt, x + (1/2) * k1x - (3/2) * k3x - (2) * k4x,
207                                           y + (1/2) * k1y - (3/2) * k3y - (2) * k4y)
208        k5x = dt * k5x
209        k5y = dt * k5y
210        --
211        if adaptive then
212            local teps = sqrt(((1/10-1/6) * k1x + (3/10) * k3x + (2/5-2/3) * k4x + (1/5 -1/6) * k5x)^2 +
213                              ((1/10-1/6) * k1y + (3/10) * k3y + (2/5-2/3) * k4y + (1/5 -1/6) * k5y)^2 )
214            local step = eps/teps
215            if teps < eps then
216                step = step^(1/4)
217                dt = 0.9 * dt * step
218            else
219                step = step^(1/3)
220                dt = 0.9 * dt * step
221                goto again
222            end
223        end
224      ::append::
225        t = t + dt
226        x = x + (1/10) * k1x + (3/10) * k3x + (2/5) * k4x + (1/5) * k5x
227        y = y + (1/10) * k1y + (3/10) * k3y + (2/5) * k4y + (1/5) * k5y
228        r = r + 1
229        if kind == "xy" then
230            result[r] = { x, y }
231        elseif kind == "tx" then
232            result[r] = { t, x }
233        else
234            result[r] = { t, y }
235        end
236        if r >= rmax then
237         -- report("pathmax is set to %i, quiting",rmax)
238            break
239        end
240      ::again::
241    end
242    if name and hash then
243        setcache(hash,name,result)
244    end
245    return result
246end
247
248mathfunctions.rungekutta = rungekutta
249
250mp.registerfunction {
251    name   = "rungekutta",
252    action = function(specification)
253        local result = rungekutta(specification)
254        if result then
255            injectpath(result)
256        else
257            injectpath { { 0, 0 } }
258        end
259    end
260}
261
Source Browser ?
