[16620] | 1 | package symexpr
|
---|
| 2 |
|
---|
| 3 | import (
|
---|
| 4 | "fmt"
|
---|
| 5 | "math"
|
---|
| 6 | )
|
---|
| 7 |
|
---|
| 8 | func (*Time) Eval(t float64, x, c, s []float64) float64 { return t }
|
---|
| 9 |
|
---|
| 10 | func (v *Var) Eval(t float64, x, c, s []float64) float64 { return x[v.P] }
|
---|
| 11 |
|
---|
| 12 | func (cnst *Constant) Eval(t float64, x, c, s []float64) float64 { return c[cnst.P] }
|
---|
| 13 |
|
---|
| 14 | func (cnst *ConstantF) Eval(t float64, x, c, s []float64) float64 { return cnst.F }
|
---|
| 15 |
|
---|
| 16 | func (sys *System) Eval(t float64, x, c, s []float64) float64 { return s[sys.P] }
|
---|
| 17 |
|
---|
| 18 | func (u *Neg) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 19 | return -1. * u.C.Eval(t, x, c, s)
|
---|
| 20 | }
|
---|
| 21 |
|
---|
| 22 | func (u *Abs) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 23 | return math.Abs(u.C.Eval(t, x, c, s))
|
---|
| 24 | }
|
---|
| 25 |
|
---|
| 26 | func (u *Sqrt) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 27 | return math.Sqrt(u.C.Eval(t, x, c, s))
|
---|
| 28 | }
|
---|
| 29 |
|
---|
| 30 | func (u *Sin) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 31 | return math.Sin(u.C.Eval(t, x, c, s))
|
---|
| 32 | }
|
---|
| 33 |
|
---|
| 34 | func (u *Cos) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 35 | return math.Cos(u.C.Eval(t, x, c, s))
|
---|
| 36 | }
|
---|
| 37 |
|
---|
| 38 | func (u *Tan) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 39 | return math.Tan(u.C.Eval(t, x, c, s))
|
---|
| 40 | }
|
---|
| 41 |
|
---|
| 42 | func (u *Exp) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 43 | return math.Exp(u.C.Eval(t, x, c, s))
|
---|
| 44 | }
|
---|
| 45 |
|
---|
| 46 | func (u *Log) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 47 | return math.Log(u.C.Eval(t, x, c, s))
|
---|
| 48 | }
|
---|
| 49 |
|
---|
| 50 | func (u *PowI) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 51 | return math.Pow(u.Base.Eval(t, x, c, s), float64(u.Power))
|
---|
| 52 | }
|
---|
| 53 |
|
---|
| 54 | func (u *PowF) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 55 | return math.Pow(u.Base.Eval(t, x, c, s), u.Power)
|
---|
| 56 | }
|
---|
| 57 |
|
---|
| 58 | func (n *PowE) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 59 | return math.Pow(n.Base.Eval(t, x, c, s), n.Power.Eval(t, x, c, s))
|
---|
| 60 | }
|
---|
| 61 |
|
---|
| 62 | func (n *Div) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 63 | return n.Numer.Eval(t, x, c, s) / n.Denom.Eval(t, x, c, s)
|
---|
| 64 | }
|
---|
| 65 |
|
---|
| 66 | func (n *Add) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 67 | ret := 0.0
|
---|
| 68 | for _, C := range n.CS {
|
---|
| 69 | if C == nil {
|
---|
| 70 | continue
|
---|
| 71 | }
|
---|
| 72 | ret += C.Eval(t, x, c, s)
|
---|
| 73 | }
|
---|
| 74 | return ret
|
---|
| 75 | }
|
---|
| 76 |
|
---|
| 77 | func (n *Mul) Eval(t float64, x, c, s []float64) float64 {
|
---|
| 78 | ret := 1.0
|
---|
| 79 | for _, C := range n.CS {
|
---|
| 80 | if C == nil {
|
---|
| 81 | continue
|
---|
| 82 | }
|
---|
| 83 | ret *= C.Eval(t, x, c, s)
|
---|
| 84 | }
|
---|
| 85 | return ret
|
---|
| 86 | }
|
---|
| 87 |
|
---|
| 88 | func RK4(eqn []Expr, ti, tj float64, x_in, x_tmp, c, s []float64) (x_out []float64) {
|
---|
| 89 | var k [32][4]float64
|
---|
| 90 | L := len(x_in)
|
---|
| 91 | h := tj - ti
|
---|
| 92 | for i := 0; i < L; i++ {
|
---|
| 93 | k[i][0] = eqn[i].Eval(ti, x_in, c, s)
|
---|
| 94 | }
|
---|
| 95 | for i := 0; i < L; i++ {
|
---|
| 96 | x_tmp[i] = x_in[i] + (h * k[i][0] / 2.0)
|
---|
| 97 | }
|
---|
| 98 | for i := 0; i < L; i++ {
|
---|
| 99 | k[i][1] = eqn[i].Eval(ti, x_tmp, c, s)
|
---|
| 100 | }
|
---|
| 101 | for i := 0; i < L; i++ {
|
---|
| 102 | x_tmp[i] = x_in[i] + (h * k[i][1] / 2.0)
|
---|
| 103 | }
|
---|
| 104 | for i := 0; i < L; i++ {
|
---|
| 105 | k[i][2] = eqn[i].Eval(ti, x_tmp, c, s)
|
---|
| 106 | }
|
---|
| 107 | for i := 0; i < L; i++ {
|
---|
| 108 | x_tmp[i] = x_in[i] + (h * k[i][2])
|
---|
| 109 | }
|
---|
| 110 | for i := 0; i < L; i++ {
|
---|
| 111 | k[i][3] = eqn[i].Eval(ti, x_tmp, c, s)
|
---|
| 112 | }
|
---|
| 113 | for i := 0; i < L; i++ {
|
---|
| 114 | x_out[i] = ((k[i][0] + 2.0*k[i][1] + 2.0*k[i][2] + k[i][3]) * (h / 6.0))
|
---|
| 115 | }
|
---|
| 116 |
|
---|
| 117 | return
|
---|
| 118 | }
|
---|
| 119 |
|
---|
| 120 | func PRK4(xn int, eqn Expr, ti, tj float64, x_in, x_out, x_tmp, c, s []float64) float64 {
|
---|
| 121 | var k [4]float64
|
---|
| 122 | L := len(x_in)
|
---|
| 123 | h := tj - ti
|
---|
| 124 | for i := 0; i < L; i++ {
|
---|
| 125 | mid := (0.5 * (x_out[i] - x_in[i]))
|
---|
| 126 | x_tmp[i] = x_in[i] + mid
|
---|
| 127 | }
|
---|
| 128 | k[0] = eqn.Eval(ti, x_in, c, s)
|
---|
| 129 | x_tmp[xn] = x_in[xn] + (h * k[0] / 2.0)
|
---|
| 130 | k[1] = eqn.Eval(ti, x_tmp, c, s)
|
---|
| 131 | x_tmp[xn] = x_in[xn] + (h * k[1] / 2.0)
|
---|
| 132 | k[2] = eqn.Eval(ti, x_tmp, c, s)
|
---|
| 133 | x_tmp[xn] = x_in[xn] + (h * k[2])
|
---|
| 134 | k[3] = eqn.Eval(ti, x_tmp, c, s)
|
---|
| 135 | return ((k[0] + 2.0*k[1] + 2.0*k[2] + k[3]) * (h / 6.0))
|
---|
| 136 | }
|
---|
| 137 |
|
---|
| 138 | func PrintPRK4(xn int, eqn Expr, ti, to float64, x_in, x_out, x_tmp, c, s []float64) float64 {
|
---|
| 139 | var k [4]float64
|
---|
| 140 | L := len(x_in)
|
---|
| 141 | h := to - ti
|
---|
| 142 | for i := 0; i < L; i++ {
|
---|
| 143 | x_tmp[i] = x_in[i] + (0.5 * (x_out[i] - x_in[i]))
|
---|
| 144 | }
|
---|
| 145 | fmt.Printf("in: %v\n", x_in)
|
---|
| 146 | fmt.Printf("out: %v\n", x_out)
|
---|
| 147 |
|
---|
| 148 | fmt.Printf("tmp: %v\n", x_tmp)
|
---|
| 149 | k[0] = eqn.Eval(ti, x_in, c, s)
|
---|
| 150 | x_tmp[xn] = x_in[xn] + (h * k[0] / 2.0)
|
---|
| 151 | fmt.Printf("tmp: %v\n", x_tmp)
|
---|
| 152 | k[1] = eqn.Eval(ti, x_tmp, c, s)
|
---|
| 153 | x_tmp[xn] = x_in[xn] + (h * k[1] / 2.0)
|
---|
| 154 | fmt.Printf("tmp: %v\n", x_tmp)
|
---|
| 155 | k[2] = eqn.Eval(ti, x_tmp, c, s)
|
---|
| 156 | x_tmp[xn] = x_in[xn] + (h * k[2])
|
---|
| 157 | fmt.Printf("tmp: %v\n", x_tmp)
|
---|
| 158 | k[3] = eqn.Eval(ti, x_tmp, c, s)
|
---|
| 159 | fmt.Printf("k: %v\n", k)
|
---|
| 160 | ans := ((k[0] + 2.0*k[1] + 2.0*k[2] + k[3]) * (h / 6.0))
|
---|
| 161 | fmt.Printf("ans: %.4f => %.4f\n\n", ans, x_out[xn]-x_in[xn])
|
---|
| 162 | return ans
|
---|
| 163 | }
|
---|