30 lines
901 B
Text
30 lines
901 B
Text
local fmt = require "fmt"
|
|
|
|
local function runge_kutta4(t0, tz, dt, y, yd)
|
|
local tn = t0
|
|
local yn = y(tn)
|
|
local z = (tz - t0) // dt
|
|
for i = 0, z do
|
|
if i % 10 == 0 then
|
|
local exact = y(tn)
|
|
local err = yn - exact
|
|
fmt.print("%4.1f %10f %10f %9f", tn, yn, exact, err)
|
|
end
|
|
if i == z then break end
|
|
local dy1 = dt * yd(tn, yn)
|
|
local dy2 = dt * yd(tn + 0.5 * dt, yn + 0.5 * dy1)
|
|
local dy3 = dt * yd(tn + 0.5 * dt, yn + 0.5 * dy2)
|
|
local dy4 = dt * yd(tn + dt, yn + dy3)
|
|
yn += (dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4) / 6.0
|
|
tn += dt
|
|
end
|
|
end
|
|
|
|
print(" T RK4 Exact Error")
|
|
print("---- --------- ---------- ---------")
|
|
local function y(t)
|
|
local x = t * t + 4.0
|
|
return x * x / 16.0
|
|
end
|
|
local yd = |t, yt| -> t * math.sqrt(yt)
|
|
runge_kutta4(0, 10, 0.1, y, yd)
|