Runge-Kutta method: Difference between revisions

import math, strformat

proc fn(t, y: float): float =

result = t * math.sqrt(y)

proc solution(t: float): float =

result = (t^2 + 4)^2 / 16

proc rk(start, stop, step: float) =

let nsteps = int(round((stop - start) / step)) + 1

let delta = (stop - start) / float(nsteps - 1)

var cur_y = 1.0

for i in 0..<nsteps:

let cur_t = start + delta * float(i)

if abs(cur_t - math.round(cur_t)) < 1e-5:

echo &"y({cur_t}) = {cur_y}, error = {solution(cur_t) - cur_y}"

let dy1 = step * fn(cur_t, cur_y)

let dy2 = step * fn(cur_t + 0.5 * step, cur_y + 0.5 * dy1)

let dy3 = step * fn(cur_t + 0.5 * step, cur_y + 0.5 * dy2)

let dy4 = step * fn(cur_t + step, cur_y + dy3)

cur_y += (dy1 + 2 * (dy2 + dy3) + dy4) / 6

rk(start = 0, stop = 10, step = 0.1)

cur_y += (dy1 + 2.0 * (dy2 + dy3) + dy4) </lang>

y(1.0) = 1.562499854278108, error = 1.457218921085968e-07

y(2.0) = 3.9999990805208, error = 9.194792003341945e-07

y(3.0) = 10.56249709043755, error = 2.909562448749625e-06

y(4.0) = 24.99999376509064, error = 6.234909363911356e-06

y(5.0) = 52.56248918030258, error = 1.081969741534294e-05

y(6.0) = 99.99998340540358, error = 1.659459641700778e-05

y(7.0) = 175.5624764822713, error = 2.351772874931157e-05

y(8.0) = 288.9999684347986, error = 3.156520142510999e-05

y(9.0) = 451.5624592768397, error = 4.07231603389846e-05

y(10.0) = 675.9999490167097, error = 5.098329029351589e-05</pre>