Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
96131 views
1
# -*- coding: utf-8 -*-
2
3
"""Copyright 2015 Roger R Labbe Jr.
4
5
FilterPy library.
6
http://github.com/rlabbe/filterpy
7
8
Documentation at:
9
https://filterpy.readthedocs.org
10
11
Supporting book at:
12
https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python
13
14
This is licensed under an MIT license. See the readme.MD file
15
for more information.
16
"""
17
18
#
19
from numpy import array, asarray, isscalar, eye, dot
20
from functools import reduce
21
22
23
def dot3(A,B,C):
24
""" Returns the matrix multiplication of A*B*C"""
25
return dot(A, dot(B,C))
26
27
def dot4(A,B,C,D):
28
""" Returns the matrix multiplication of A*B*C*D"""
29
return dot(A, dot(B, dot(C,D)))
30
31
32
def dotn(*args):
33
""" Returns the matrix multiplication of 2 or more matrices"""
34
return reduce(dot, args)
35
36
37
def runge_kutta4(y, x, dx, f):
38
"""computes 4th order Runge-Kutta for dy/dx.
39
40
**Parameters**
41
42
y : scalar
43
Initial/current value for y
44
x : scalar
45
Initial/current value for x
46
dx : scalar
47
difference in x (e.g. the time step)
48
f : ufunc(y,x)
49
Callable function (y, x) that you supply to compute dy/dx for
50
the specified values.
51
52
"""
53
54
k1 = dx * f(y, x)
55
k2 = dx * f(y + 0.5*k1, x + 0.5*dx)
56
k3 = dx * f(y + 0.5*k2, x + 0.5*dx)
57
k4 = dx * f(y + k3, x + dx)
58
59
return y + (k1 + 2*k2 + 2*k3 + k4) / 6.
60
61
62
def setter(value, dim_x, dim_y):
63
""" Returns a copy of 'value' as an numpy.array with dtype=float. Throws
64
exception if the array is not dimensioned correctly. Value may be any
65
type which converts to numpy.array (list, np.array, np.matrix, etc)
66
"""
67
v = array(value, dtype=float)
68
if v.shape != (dim_x, dim_y):
69
raise Exception('must have shape ({},{})'.format(dim_x, dim_y))
70
return v
71
72
def setter_1d(value, dim_x):
73
""" Returns a copy of 'value' as an numpy.array with dtype=float. Throws
74
exception if the array is not dimensioned correctly. Value may be any
75
type which converts to numpy.array (list, np.array, np.matrix, etc)
76
"""
77
78
v = array(value, dtype=float)
79
shape = v.shape
80
if shape[0] != (dim_x) or v.ndim > 2 or (v.ndim==2 and shape[1] != 1):
81
raise Exception('has shape {}, must have shape ({},{})'.format(shape, dim_x, 1))
82
return v
83
84
85
def setter_scalar(value, dim_x):
86
""" Returns a copy of 'value' as an numpy.array with dtype=float. Throws
87
exception if the array is not dimensioned correctly. Value may be any
88
type which converts to numpy.array (list, np.array, np.matrix, etc),
89
or a scalar, in which case we create a diagonal matrix with each
90
diagonal element == value.
91
"""
92
if isscalar(value):
93
v = eye(dim_x) * value
94
else:
95
v = asarray(value, dtype=float)
96
97
if v is value:
98
v = value.copy()
99
if v.shape != (dim_x, dim_x):
100
raise Exception('must have shape ({},{})'.format(dim_x, dim_x))
101
return v
102
103