/
fun_orthopoly01.go
107 lines (95 loc) · 3.03 KB
/
fun_orthopoly01.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
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
// 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/fun"
"github.com/cpmech/gosl/io"
"github.com/cpmech/gosl/plt"
"github.com/cpmech/gosl/utl"
)
func main() {
// Reference:
// [1] Abramowitz M, Stegun IA (1972) Handbook of Mathematical Functions with Formulas,
// Graphs, and Mathematical Tables. U.S. Department of Commerce, NIST
// number of points
npts := 501
// Jacobi polynomials Fig 22.1
N, α, β := 5, 1.5, -0.5
jacobi := fun.NewGeneralOrthoPoly(fun.OpolyJacobiKind, N, α, β)
plt.Reset(true, &plt.A{Prop: 1.5})
X := utl.LinSpace(-1, 1, npts)
Y := make([]float64, len(X))
for n := 0; n <= 5; n++ {
for i := 0; i < len(X); i++ {
Y[i] = jacobi.P(n, X[i])
}
plt.Plot(X, Y, &plt.A{L: io.Sf("$P_{%d}^{(%g,%g)}$", n, α, β), NoClip: true})
}
plt.Cross(0, 0, nil)
plt.Equal()
plt.AxisYrange(-1, 3.3)
plt.HideAllBorders()
plt.Gll("$x$", io.Sf("$P_n^{(%g,%g)}$", α, β), &plt.A{LegOut: true, LegNcol: 3})
plt.Save("/tmp/gosl", "as-fig-22-1")
// Chebyshev1 polynomials Fig 22.6
chebyshev1 := fun.NewGeneralOrthoPoly(fun.OpolyCheby1Kind, N, 0, 0)
plt.Reset(true, &plt.A{Prop: 0.8})
for n := 1; n <= 5; n++ {
for i := 0; i < len(X); i++ {
Y[i] = chebyshev1.P(n, X[i])
}
plt.Plot(X, Y, &plt.A{L: io.Sf("$P_{%d}$", n), NoClip: true})
}
plt.Cross(0, 0, nil)
plt.HideAllBorders()
plt.Gll("$x$", "$P_n$", &plt.A{LegOut: true, LegNcol: 3})
plt.Save("/tmp/gosl", "as-fig-22-6")
// Chebyshev2 polynomials Fig 22.7
chebyshev2 := fun.NewGeneralOrthoPoly(fun.OpolyCheby2Kind, N, 0, 0)
plt.Reset(true, &plt.A{Prop: 0.8})
for n := 1; n <= 5; n++ {
for i := 0; i < len(X); i++ {
Y[i] = chebyshev2.P(n, X[i])
}
plt.Plot(X, Y, &plt.A{L: io.Sf("$P_{%d}$", n), NoClip: true})
}
plt.Cross(0, 0, nil)
plt.HideAllBorders()
plt.AxisYrange(-3, 5.5)
plt.Gll("$x$", "$P_n$", &plt.A{LegOut: true, LegNcol: 3})
plt.Save("/tmp/gosl", "as-fig-22-7")
// Legendre polynomials Fig 22.8
legendre := fun.NewGeneralOrthoPoly(fun.OpolyLegendreKind, N, 0, 0)
plt.Reset(true, &plt.A{Prop: 1.0})
for n := 0; n <= 5; n++ {
for i := 0; i < len(X); i++ {
Y[i] = legendre.P(n, X[i])
}
plt.Plot(X, Y, &plt.A{L: io.Sf("$P_{%d}$", n), NoClip: true})
}
plt.Cross(0, 0, nil)
plt.Equal()
plt.AxisYrange(-0.5, 1.0)
plt.HideAllBorders()
plt.Gll("$x$", "$P_n$", &plt.A{LegOut: true, LegNcol: 3})
plt.Save("/tmp/gosl", "as-fig-22-8")
// Hermite polynomials Fig 22.10
hermite := fun.NewGeneralOrthoPoly(fun.OpolyHermiteKind, N, 0, 0)
plt.Reset(true, &plt.A{Prop: 0.8})
X = utl.LinSpace(0, 4, npts)
for n := 2; n <= 5; n++ {
den := math.Pow(float64(n), 3)
for i := 0; i < len(X); i++ {
Y[i] = hermite.P(n, X[i]) / den
}
plt.Plot(X, Y, &plt.A{L: io.Sf("$P_{%d}/%g$", n, den), NoClip: true})
}
plt.Cross(0, 0, nil)
plt.AxisYrange(-1, 8.5)
plt.HideAllBorders()
plt.Gll("$x$", "$P_n$", &plt.A{LegOut: true, LegNcol: 3})
plt.Save("/tmp/gosl", "as-fig-22-10")
}