forked from vlang/vsl
-
Notifications
You must be signed in to change notification settings - Fork 0
/
diff.v
122 lines (118 loc) · 3.33 KB
/
diff.v
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
module diff
import vsl.func
import vsl.internal.prec
import math
pub fn backward(f func.Fn, x f64) (f64, f64) {
/*
Construct a divided difference table with a fairly large step
* size to get a very rough estimate of f''. Use this to estimate
* the step size which will minimize the error in calculating f'.
*/
mut h := prec.sqrt_f64_epsilon
mut a := []f64{}
mut d := []f64{}
mut k := 0
mut i := 0
/*
Algorithm based on description on pg. 204 of Conte and de Boor
* (CdB) - coefficients of Newton form of polynomial of degree 2.
*/
for i = 0; i < 3; i++ {
a << x + (f64(i) - 2.0) * h
d << f.eval(a[i])
}
for k = 1; k < 4; k++ {
for i = 0; i < 3 - k; i++ {
d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i])
}
}
/*
Adapt procedure described on pg. 282 of CdB to find best value of
* step size.
*/
mut a2 := math.abs(d[0] + d[1] + d[2])
if a2 < 100.0 * prec.sqrt_f64_epsilon {
a2 = 100.0 * prec.sqrt_f64_epsilon
}
h = math.sqrt(prec.sqrt_f64_epsilon / (2.0 * a2))
if h > 100.0 * prec.sqrt_f64_epsilon {
h = 100.0 * prec.sqrt_f64_epsilon
}
return (f.eval(x) - f.eval(x - h)) / h, math.abs(10.0 * a2 * h)
}
pub fn forward(f func.Fn, x f64) (f64, f64) {
/*
Construct a divided difference table with a fairly large step
* size to get a very rough estimate of f''. Use this to estimate
* the step size which will minimize the error in calculating f'.
*/
mut h := prec.sqrt_f64_epsilon
mut a := []f64{}
mut d := []f64{}
mut k := 0
mut i := 0
/*
Algorithm based on description on pg. 204 of Conte and de Boor
* (CdB) - coefficients of Newton form of polynomial of degree 2.
*/
for i = 0; i < 3; i++ {
a << x + f64(i) * h
d << f.eval(a[i])
}
for k = 1; k < 4; k++ {
for i = 0; i < 3 - k; i++ {
d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i])
}
}
/*
Adapt procedure described on pg. 282 of CdB to find best value of
* step size.
*/
mut a2 := math.abs(d[0] + d[1] + d[2])
if a2 < 100.0 * prec.sqrt_f64_epsilon {
a2 = 100.0 * prec.sqrt_f64_epsilon
}
h = math.sqrt(prec.sqrt_f64_epsilon / (2.0 * a2))
if h > 100.0 * prec.sqrt_f64_epsilon {
h = 100.0 * prec.sqrt_f64_epsilon
}
return (f.eval(x + h) - f.eval(x)) / h, math.abs(10.0 * a2 * h)
}
pub fn central(f func.Fn, x f64) (f64, f64) {
/*
Construct a divided difference table with a fairly large step
* size to get a very rough estimate of f'''. Use this to estimate
* the step size which will minimize the error in calculating f'.
*/
mut h := prec.sqrt_f64_epsilon
mut a := []f64{}
mut d := []f64{}
mut k := 0
mut i := 0
/*
Algorithm based on description on pg. 204 of Conte and de Boor
* (CdB) - coefficients of Newton form of polynomial of degree 3.
*/
for i = 0; i < 4; i++ {
a << x + (f64(i) - 2.0) * h
d << f.eval(a[i])
}
for k = 1; k < 5; k++ {
for i = 0; i < 3 - k; i++ {
d[i] = (d[i + 1] - d[i]) / (a[i + k] - a[i])
}
}
/*
Adapt procedure described on pg. 282 of CdB to find best value of
* step size.
*/
mut a3 := math.abs(d[0] + d[1] + d[2] + d[3])
if a3 < 100.0 * prec.sqrt_f64_epsilon {
a3 = 100.0 * prec.sqrt_f64_epsilon
}
h = math.pow(prec.sqrt_f64_epsilon / (2.0 * a3), 1.0 / 3.0)
if h > 100.0 * prec.sqrt_f64_epsilon {
h = 100.0 * prec.sqrt_f64_epsilon
}
return (f.eval(x + h) - f.eval(x - h)) / (2.0 * h), math.abs(100.0 * a3 * h * h)
}