1
2
3
4 from __future__ import division, absolute_import
5
6 __author__ = "Eduardo dos Santos Pereira"
7 __email__ = "pereira.somoza@gmail.com"
8 __credits__ = ["Eduardo dos Santos Pereira"]
9 __license__ = "GPLV3"
10 __version__ = "1.0.1"
11 __maintainer__ = "Eduardo dos Santos Pereira"
12 __status__ = "Stable"
13
14
15 """4th-order Runge-Kutta method for solving the initial value problem
16 X,Y = integrate(F,x,y,xStop,h).
17 4th-order Runge-Kutta method for solving the
18 initial value problem { y} ’ = { F(x,{ y} )} , where
19 { y} = { y[0],y[1],...y[n-1]} .
20 x,y = initial conditions.
21 xStop = terminal value of x.
22 h = increment of x used in integration.
23 F = user-supplied function that returns the
24 array F(x,y) = { y’[0],y’[1],...,y’[n-1]} .
25
26
27 This file is part of pystar.
28 copyright : Eduardo dos Santos Pereira
29
30 pystar is free software: you can redistribute it and/or modify
31 it under the terms of the GNU General Public License as published by
32 the Free Software Foundation, either version 3 of the License.
33 pystar is distributed in the hope that it will be useful,
34 but WITHOUT ANY WARRANTY; without even the implied warranty of
35 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 GNU General Public License for more details.
37
38 You should have received a copy of the GNU General Public License
39 along with Foobar. If not, see <http://www.gnu.org/licenses/>.
40
41 """
42
43 from numpy import array
44
45
47
48 def run_kut4(F, x, y, h):
49
50 K0 = h * F(x, y)
51 K1 = h * F(x + h / 2.0, y + K0 / 2.0)
52 K2 = h * F(x + h / 2.0, y + K1 / 2.0)
53 K3 = h * F(x + h, y + K2)
54
55 return (K0 + 2.0 * K1 + 2.0 * K2 + K3) / 6.0
56
57 X = []
58 Y = []
59 X.append(x)
60 Y.append(y)
61
62 while x < xStop:
63 h = min(h, xStop - x)
64 y = y + run_kut4(F, x, y, h)
65 x = x + h
66
67 X.append(x)
68 Y.append(y)
69
70 return array(X), array(Y)
71