3
import "github.com/ready-steady/linear/system"
15
func NewCubic(x, y []float64) *Cubic {
19
dx := make([]float64, nn-1)
20
dydx := make([]float64, (nn-1)*nd)
21
for i := 0; i < (nn - 1); i++ {
23
for j := 0; j < nd; j++ {
24
dydx[i*nd+j] = (y[(i+1)*nd+j] - y[i*nd+j]) / dx[i]
32
s.Nodes = []float64{x[0], x[1]}
33
s.Weights = make([]float64, nd*4)
34
for j := 0; j < nd; j++ {
35
s.Weights[j*4+2] = dydx[j]
36
s.Weights[j*4+3] = y[j]
39
s.Nodes = []float64{x[0], x[2]}
40
s.Weights = make([]float64, nd*4)
41
for j := 0; j < nd; j++ {
42
c1 := (dydx[nd+j] - dydx[j]) / (x[2] - x[0])
44
s.Weights[j*4+2] = dydx[j] - c1*dx[0]
45
s.Weights[j*4+3] = y[j]
49
xe := x[nn-1] - x[nn-3]
51
a := make([]float64, nn)
52
for i := 0; i < (nn - 2); i++ {
57
b := make([]float64, nn)
59
for i := 1; i < (nn - 1); i++ {
60
b[i] = 2 * (dx[i] + dx[i-1])
64
c := make([]float64, nn)
66
for i := 2; i < nn; i++ {
70
d := make([]float64, nd*nn)
71
for j := 0; j < nd; j++ {
72
d[j*nn] = ((dx[0]+2*xb)*dx[1]*dydx[j] + dx[0]*dx[0]*dydx[nd+j]) / xb
73
for i := 1; i < (nn - 1); i++ {
74
d[j*nn+i] = 3 * (dx[i]*dydx[(i-1)*nd+j] + dx[i-1]*dydx[i*nd+j])
76
d[j*nn+nn-1] = (dx[nn-2]*dx[nn-2]*dydx[(nn-3)*nd+j] +
77
(2*xe+dx[nn-2])*dx[nn-3]*dydx[(nn-2)*nd+j]) / xe
80
slopes := system.ComputeTridiagonal(a, b, c, d)
82
s.Nodes = make([]float64, nn)
85
s.Weights = make([]float64, (nn-1)*nd*4)
86
for i, k := 0, 0; i < (nn - 1); i++ {
87
for j := 0; j < nd; j++ {
88
α := (dydx[i*nd+j] - slopes[j*nn+i]) / dx[i]
89
β := (slopes[j*nn+i+1] - dydx[i*nd+j]) / dx[i]
90
s.Weights[k] = (β - α) / dx[i]
91
s.Weights[k+1] = 2*α - β
92
s.Weights[k+2] = slopes[j*nn+i]
93
s.Weights[k+3] = y[i*nd+j]
105
func (s *Cubic) Evaluate(x []float64) []float64 {
106
nn, np := len(s.Nodes), len(x)
107
nd := len(s.Weights) / (4 * (nn - 1))
109
y := make([]float64, np*nd)
111
for i, l := 0, 0; i < np; i++ {
112
for x[i] > s.Nodes[l+1] && l < (nn-2) {
116
z := x[i] - s.Nodes[l]
117
for j, k := 0, l*nd*4; j < nd; j++ {
118
y[i*nd+j] = z*(z*(z*s.Weights[k]+s.Weights[k+1])+s.Weights[k+2]) + s.Weights[k+3]