41 lines
1.2 KiB
Text
41 lines
1.2 KiB
Text
function rk4(f, t0, y0, t1, n) {
|
|
h = (t1-t0)/(n-1)
|
|
a = J(n, 2, 0)
|
|
a[1, 1] = t = t0
|
|
a[1, 2] = y = y0
|
|
for (i=2; i<=n; i++) {
|
|
k1 = h*(*f)(t, y)
|
|
k2 = h*(*f)(t+0.5*h, y+0.5*k1)
|
|
k3 = h*(*f)(t+0.5*h, y+0.5*k2)
|
|
k4 = h*(*f)(t+h, y+k3)
|
|
t = t+h
|
|
y = y+(k1+2*k2+2*k3+k4)/6
|
|
a[i, 1] = t
|
|
a[i, 2] = y
|
|
}
|
|
return(a)
|
|
}
|
|
|
|
function f(t, y) {
|
|
return(t*sqrt(y))
|
|
}
|
|
|
|
a = rk4(&f(), 0, 1, 10, 101)
|
|
t = a[., 1]
|
|
a = a, a[., 2]:-(t:^2:+4):^2:/16
|
|
a[range(1,101,10), .]
|
|
|
|
1 2 3
|
|
+----------------------------------------------+
|
|
1 | 0 1 0 |
|
|
2 | 1 1.562499854 -1.45722e-07 |
|
|
3 | 2 3.999999081 -9.19479e-07 |
|
|
4 | 3 10.56249709 -2.90956e-06 |
|
|
5 | 4 24.99999377 -6.23491e-06 |
|
|
6 | 5 52.56248918 -.0000108197 |
|
|
7 | 6 99.99998341 -.0000165946 |
|
|
8 | 7 175.5624765 -.0000235177 |
|
|
9 | 8 288.9999684 -.0000315652 |
|
|
10 | 9 451.5624593 -.0000407232 |
|
|
11 | 10 675.999949 -.0000509833 |
|
|
+----------------------------------------------+
|