-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdemo_pspline_1d.py
127 lines (98 loc) · 2.85 KB
/
demo_pspline_1d.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
#!/usr/bin/env python
"""
1-D spline in float64 precision
"""
import sys, time
import numpy as np
import matplotlib.pyplot as plt
from pypspline3.pspline_1d import pspline
eps = 1.e-6
n1 = 11
bcs1 = (0,0)
x1min, x1max = 0., 1.
x1 = np.arange(x1min, x1max+eps, (x1max-x1min)/float(n1-1))
tic = time.time()
#########
f = x1**3
#########
toc = time.time()
print("function evaluations: time->%10.1f secs" % (toc-tic))
tic = time.time()
spl = pspline(x1)
# may set BCs if not-a-knot
spl.setup(f)
toc = time.time()
print("init/setup: %d original grid nodes time->%10.1f secs" % \
(n1, toc-tic))
# save/load is not considered here
# new mesh
n1 = 2*n1-1
x2 = np.arange(x1min, x1max+eps, (x1max-x1min)/float(n1-1))
fexact = x2**3
# point interpolation
nint = n1
# all_fi = np.zeros(n1)
error = 0
all_fi = np.zeros(n1)
tic = time.time()
for i1 in range(n1):
fi, ier, iwarn = spl.interp_point(x2[i1])
error += (fi - fexact[i1])**2
all_fi[i1] = fi
toc = time.time()
error /= nint
error = np.sqrt(error)
print("interp_point: %d evaluations (error=%g) ier=%d iwarn=%d time->%10.1f secs" % \
(nint, error, ier, iwarn, toc-tic))
plt.figure()
plt.plot(x1, f, "o-", label="orig")
plt.plot(x2, all_fi, ".--", label="interp")
plt.legend(loc="upper left")
plt.show()
# array interpolation
tic = time.time()
fi, ier, iwarn = spl.interp_array(x2)
toc = time.time()
error = np.sum(np.sum(np.sum((fi-fexact)**2)))/nint
print("interp_array: %d evaluations (error=%g) ier=%d iwarn=%d time->%10.1f secs" % \
(nint, error, ier, iwarn, toc-tic))
## df/dx
fexact = 3*x2**2
# point df/dx
tic = time.time()
error = 0
for i1 in range(n1):
fi, ier, iwarn = spl.derivative_point(1, x2[i1])
error += (fi - fexact[i1])**2
toc = time.time()
error /= nint
error = np.sqrt(error)
print("derivative_point df/dx: %d evaluations (error=%g) ier=%d iwarn=%d time->%10.1f secs" % \
(nint, error, ier, iwarn, toc-tic))
# array df/dx
tic = time.time()
fi, ier, iwarn = spl.derivative_array(1, x2)
toc = time.time()
error = np.sum(np.sum(np.sum((fi-fexact)**2)))/nint
print("derivative_array df/dx: %d evaluations (error=%g) ier=%d iwarn=%d time->%10.1f secs" % \
(nint, error, ier, iwarn, toc-tic))
## d^2f/dx^2
fexact = 3*2*x2
# point d^2f/dx^2
tic = time.time()
error = 0
for i1 in range(n1):
fi, ier, iwarn = spl.derivative_point(2, x2[i1])
error += (fi - fexact[i1])**2
toc = time.time()
error /= nint
error = np.sqrt(error)
print("derivative_point d^2f/dx^2: %d evaluations (error=%g) ier=%d iwarn=%d time->%10.1f secs" % \
(nint, error, ier, iwarn, toc-tic))
# array d^2f/dx^2
tic = time.time()
fi, ier, iwarn = spl.derivative_array(2, x2)
toc = time.time()
error = np.sum(np.sum(np.sum((fi-fexact)**2)))/nint
print("derivative_array d^2f/dx^2: %d evaluations (error=%g) ier=%d iwarn=%d time->%10.1f secs" % \
(nint, error, ier, iwarn, toc-tic))