forked from cpmech/gosl
/
num_deriv01.go
62 lines (49 loc) · 1.59 KB
/
num_deriv01.go
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
// Copyright 2016 The Gosl Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
// +build ignore
package main
import (
"math"
"github.com/cpmech/gosl/chk"
"github.com/cpmech/gosl/io"
"github.com/cpmech/gosl/num"
"github.com/cpmech/gosl/plt"
"github.com/cpmech/gosl/utl"
)
func main() {
// define function and derivative function
y_fcn := func(x float64) float64 { return math.Sin(x) }
dydx_fcn := func(x float64) float64 { return math.Cos(x) }
d2ydx2_fcn := func(x float64) float64 { return -math.Sin(x) }
// run test for 11 points
X := utl.LinSpace(0, 2*math.Pi, 11)
io.Pf(" %8s %23s %23s %23s\n", "x", "analytical", "numerical", "error")
for _, x := range X {
// analytical derivatives
dydx_ana := dydx_fcn(x)
d2ydx2_ana := d2ydx2_fcn(x)
// numerical derivative: dydx
dydx_num, _ := num.DerivCen5(x, 1e-3, func(t float64) (float64, error) {
return y_fcn(t), nil
})
// numerical derivative d2ydx2
d2ydx2_num, _ := num.DerivCen5(x, 1e-3, func(t float64) (float64, error) {
return dydx_fcn(t), nil
})
// check
chk.PrintAnaNum(io.Sf("dy/dx @ %.6f", x), 1e-10, dydx_ana, dydx_num, true)
chk.PrintAnaNum(io.Sf("d²y/dx² @ %.6f", x), 1e-10, d2ydx2_ana, d2ydx2_num, true)
}
// generate 101 points for plotting
X = utl.LinSpace(0, 2*math.Pi, 101)
Y := make([]float64, len(X))
for i, x := range X {
Y[i] = y_fcn(x)
}
// plot
plt.Reset(false, nil)
plt.Plot(X, Y, &plt.A{C: "b", M: ".", Me: 10, L: "y(x)=sin(x)"})
plt.Gll("x", "y", nil)
plt.Save("/tmp/gosl", "num_deriv01")
}