-
Notifications
You must be signed in to change notification settings - Fork 0
/
cutility.py
143 lines (116 loc) · 4.81 KB
/
cutility.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
'''
math helper functions
Franz Ludwig Kostelezky, 2021
'''
import numpy as np
from itertools import product
def finite_difference_derivate_3_point(series, _):
derivate = - np.roll(series, 1) + np.roll(series, -1)
derivate = derivate / 2
return derivate
def finite_difference_derivate_5_point(series, _):
derivate = - np.roll(series, -2) + 8 * np.roll(series, -1) - 8 * np.roll(series, 1) + np.roll(series, 2)
derivate = derivate / 12
return derivate
def finite_difference_derivate_7_point(series, _):
derivate = - np.roll(series, 3) + 9 * np.roll(series, 2) - 45 * np.roll(series, 1) + 45 * np.roll(series, -1) \
- 9 * np.roll(series, -2) + np.roll(series, -3)
derivate = derivate / 60
return derivate
def finite_difference_derivate_9_point(series, _):
derivate = - 3 * np.roll(series, 4) - 32 * np.roll(series, 3) + 168 * np.roll(series, 2) - 672 * np.roll(series, 1) \
+ 672 * np.roll(series, -1) - 168 * np.roll(series, -2) + 32 * np.roll(series, -3) - \
3 * np.roll(series, -4)
derivate = derivate / 840
return derivate
def second_order_upwind(series, dx=1):
'''Returns the 1D second order upwind derivate of a one dimensional
time series using reflecting boundary conditions.
'''
series = np.array(series)
#dx = 1
d_pos = (- 3 * series \
+ 4 * np.roll(series, shift=-1, axis=0) \
- np.roll(series, shift=-2, axis=0)
) / (2 * dx)
d_neg = (+ 3 * series \
- 4 * np.roll(series, shift=1, axis=0) \
+ np.roll(series, shift=2, axis=0)
) / (2 * dx)
derivate = d_pos
derivate[-3::] = d_neg[-3::]
return derivate
def first_order_upwind(series, dx=1):
''' first order upwind first order 1d derivate
'''
series = np.array(series)
#dx = 1
d_pos = (series - np.roll(series, shift=1, axis=0)) / dx
d_neg = (np.roll(series, shift=-1, axis=0) - series) / dx
derivate = d_pos
derivate[-2::] = d_neg[-2::]
return derivate
def third_order_upwind(series, dx=1):
''' third order upwind third order 1d derivate
'''
series = np.array(series)
#dx = 1
d_pos = (- 2 * np.roll(series, shift=1, axis=0) \
- 3 * series \
+ 6 * np.roll(series, shift=-1, axis=0) \
- np.roll(series, shift=-2, axis=0)
) / (6 * dx)
d_neg = (+ 2 * np.roll(series, shift=-1, axis=0) \
+ 3 * series \
- 6 * np.roll(series, shift=1, axis=0) \
+ np.roll(series, shift=2, axis=0)
) / (6 * dx)
derivate = d_pos
derivate[-4::] = d_neg[-4::]
return derivate
def polynominal(dimension, grade):
''' returns the exponents of a polynominal
of a given dimension to a given grade.
'''
# terminal condition
if grade == 1:
return np.identity(dimension)
# get all possible combinations of grade x dimension
tmp = product(range(grade + 1), repeat=dimension)
tmp = list(tmp)
# remove all which do not match grade
tmp_ = []
for i in range(len(tmp)):
if np.sum(tmp[i]) == grade:
tmp_.append(list(tmp[i]))
# convert to full numpy array
tmp_ = np.asarray([np.asarray(el) for el in tmp_])
return np.append(polynominal(dimension, grade - 1), tmp_.T, axis=1)
#def _polynominal(dimension, grade):
# ''' This function is obsolete. it was used to verify polynominal.
# '''
# if dimension != 3:
# return polynominal(dimension, grade)
#
# assert(dimension == 3)
# assert(grade <= 4)
# print('using 3d static - warning: obsolete')
#
# if grade == 1:
# r = [[1., 0., 0.],
# [0., 1., 0.],
# [0., 0., 1.]]
# if grade == 2: # v
# r = [[1., 0., 0., 2., 1., 1., 0., 0., 0.],
# [0., 1., 0., 0., 1., 0., 2., 0., 1.],
# [0., 0., 1., 0., 0., 1., 0., 2., 1.]]
# if grade == 3: # v
# r = [[1., 0., 0., 2., 1., 1., 0., 0., 0., 3., 2., 2., 1., 1., 1., 0., 0., 0., 0.],
# [0., 1., 0., 0., 1., 0., 2., 0., 1., 0., 1., 0., 2., 0., 1., 3., 0., 1., 2.],
# [0., 0., 1., 0., 0., 1., 0., 2., 1., 0., 0., 1., 0., 2., 1., 0., 3., 2., 1.]]
# if grade == 4: # v
# r = [[1., 0., 0., 2., 1., 1., 0., 0., 0., 3., 2., 2., 1., 1., 1., 0., 0., 0., 0., 4., 3., 3., 2., 2., 2., 1., 1., 1., 1., 0., 0., 0., 0., 0.],
# [0., 1., 0., 0., 1., 0., 2., 0., 1., 0., 1., 0., 2., 0., 1., 3., 0., 1., 2., 0., 1., 0., 2., 0., 1., 3., 0., 2., 1., 3., 1., 2., 4., 0.],
# [0., 0., 1., 0., 0., 1., 0., 2., 1., 0., 0., 1., 0., 2., 1., 0., 3., 2., 1., 0., 0., 1., 0., 2., 1., 0., 3., 1., 2., 1., 3., 2., 0., 4.]]
#
# return np.asarray(r)