23 lines
616 B
Text
23 lines
616 B
Text
F rk4(f, x0, y0, x1, n)
|
||
V vx = [0.0] * (n + 1)
|
||
V vy = [0.0] * (n + 1)
|
||
V h = (x1 - x0) / Float(n)
|
||
V x = x0
|
||
V y = y0
|
||
vx[0] = x
|
||
vy[0] = y
|
||
L(i) 1..n
|
||
V k1 = h * f(x, y)
|
||
V k2 = h * f(x + 0.5 * h, y + 0.5 * k1)
|
||
V k3 = h * f(x + 0.5 * h, y + 0.5 * k2)
|
||
V k4 = h * f(x + h, y + k3)
|
||
vx[i] = x = x0 + i * h
|
||
vy[i] = y = y + (k1 + k2 + k2 + k3 + k3 + k4) / 6
|
||
R (vx, vy)
|
||
|
||
F f(Float x, Float y) -> Float
|
||
R x * sqrt(y)
|
||
|
||
V (vx, vy) = rk4(f, 0.0, 1.0, 10.0, 100)
|
||
L(x, y) zip(vx, vy)[(0..).step(10)]
|
||
print(‘#2.1 #4.5 #2.8’.format(x, y, y - (4 + x * x) ^ 2 / 16))
|