|
| 1 | +# Reference : Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application |
| 2 | +using IntervalArithmetic, ForwardDiff |
| 3 | +import Base: +, -, *, /, ^, sqrt, exp, log, sin, cos, tan, asin, acos, atan |
| 4 | +import IntervalArithmetic: mid, interval |
| 5 | + |
| 6 | +function slope(f::Function, x::Interval, c::Number) |
| 7 | + f(slope_var(x, c)).s |
| 8 | +end |
| 9 | + |
| 10 | +struct Slope{T} |
| 11 | + x::Interval{T} # Interval on which slope is evaluated |
| 12 | + c::Interval{T} # Point about which slope is evaluated (Interval to get bounded rounding errors) |
| 13 | + s::Interval{T} # Variable which propogates the slope information |
| 14 | +end |
| 15 | + |
| 16 | +Slope(c) = Slope(c, c, 0) |
| 17 | +Slope(a, b, c) = Slope(promote(convert(Interval, a), b, c)...) |
| 18 | + |
| 19 | +function slope_var(v::Number) |
| 20 | + Slope(v, v, 1) |
| 21 | +end |
| 22 | + |
| 23 | +function slope_var(v::Interval, c::Number) |
| 24 | + Slope(v, c, 1) |
| 25 | +end |
| 26 | + |
| 27 | +function interval(u::Slope) |
| 28 | + u.x |
| 29 | +end |
| 30 | + |
| 31 | +function mid(u::Slope) |
| 32 | + u.c |
| 33 | +end |
| 34 | + |
| 35 | +function slope(u::Slope) |
| 36 | + u.s |
| 37 | +end |
| 38 | + |
| 39 | +function +(u::Slope, v::Slope) |
| 40 | + Slope(u.x + v.x, u.c + v.c, u.s + v.s) |
| 41 | +end |
| 42 | + |
| 43 | +function -(u::Slope, v::Slope) |
| 44 | + Slope(u.x - v.x, u.c - v.c, u.s - v.s) |
| 45 | +end |
| 46 | + |
| 47 | +function *(u::Slope, v::Slope) |
| 48 | + Slope(u.x * v.x, u.c * v.c, u.s * v.c + u.x * v.s) |
| 49 | +end |
| 50 | + |
| 51 | +function /(u::Slope, v::Slope) |
| 52 | + Slope(u.x / v.x, u.c / v.c, (u.s - (u.c / v.c) * v.s) / v.x) |
| 53 | +end |
| 54 | + |
| 55 | +function +(u, v::Slope) |
| 56 | + Slope(u + v.x, u + v.c, v.s) |
| 57 | +end |
| 58 | + |
| 59 | +function -(u, v::Slope) |
| 60 | + Slope(u - v.x, u - v.c, -v.s) |
| 61 | +end |
| 62 | + |
| 63 | +function *(u, v::Slope) |
| 64 | + Slope(u * v.x, u * v.c, u * v.s) |
| 65 | +end |
| 66 | + |
| 67 | +function /(u, v::Slope) |
| 68 | + Slope(u / v.x, u / v.c, -(u / v.c) * (v.s / v.x)) |
| 69 | +end |
| 70 | + |
| 71 | ++(v::Slope, u) = u + v |
| 72 | + |
| 73 | +*(v::Slope, u) = u * v |
| 74 | + |
| 75 | +function -(u::Slope, v) |
| 76 | + Slope(u.x - v, u.c - v, u.s) |
| 77 | +end |
| 78 | + |
| 79 | +function -(u::Slope) |
| 80 | + Slope(-u.x, -u.c, -u.s) |
| 81 | +end |
| 82 | + |
| 83 | +function /(u::Slope, v) |
| 84 | + Slope(u.x / v, u.c / v, u.s / v) |
| 85 | +end |
| 86 | + |
| 87 | +function sqr(u::Slope) |
| 88 | + Slope(u.x ^ 2, u.c ^ 2, (u.x + u.c) * u.s) |
| 89 | +end |
| 90 | + |
| 91 | +function ^(u::Slope, k::Integer) |
| 92 | + if k == 0 |
| 93 | + return Slope(1) |
| 94 | + elseif k == 1 |
| 95 | + return u |
| 96 | + elseif k == 2 |
| 97 | + return sqr(u) |
| 98 | + else |
| 99 | + hxi = interval(u.x.lo) ^ k |
| 100 | + hxs = interval(u.x.hi) ^ k |
| 101 | + hx = hull(hxi, hxs) |
| 102 | + |
| 103 | + if (k % 2 == 0) && (0 ∈ u.x) |
| 104 | + hx = interval(0, hx.hi) |
| 105 | + end |
| 106 | + |
| 107 | + hc = u.c ^ k |
| 108 | + |
| 109 | + i = u.x.lo - u.c.lo |
| 110 | + s = u.x.hi - u.c.hi |
| 111 | + |
| 112 | + if ((i == 0) || (s == 0) || (k % 2 == 1 && zero(u.x) ⪽ u.x)) |
| 113 | + h1 = k * (u.x ^ (k - 1)) |
| 114 | + else |
| 115 | + if k % 2 == 0 || u.x.lo >= 0 |
| 116 | + h1 = interval((hxi.hi - hc.lo) / i, (hxs.hi - hc.lo) / s) |
| 117 | + else |
| 118 | + h1 = interval((hxs.lo - hc.hi) / s, (hxi.lo - hc.hi) / i) |
| 119 | + end |
| 120 | + end |
| 121 | + return Slope(hx, hc, h1 * u.s) |
| 122 | + end |
| 123 | +end |
| 124 | + |
| 125 | +function sqrt(u::Slope) |
| 126 | + Slope(sqrt(u.x), sqrt(u.c), u.s / (sqrt(u.x) + sqrt(u.c))) |
| 127 | +end |
| 128 | + |
| 129 | +function exp(u::Slope) |
| 130 | + hx = exp(u.x) |
| 131 | + hc = exp(u.c) |
| 132 | + |
| 133 | + i = u.x.lo - u.c.lo |
| 134 | + s = u.x.hi - u.c.hi |
| 135 | + |
| 136 | + if (i == 0 || s == 0) |
| 137 | + h1 = hx |
| 138 | + else |
| 139 | + h1 = interval((hx.lo - hc.lo) / i, (hx.hi - hc.hi) / s) |
| 140 | + end |
| 141 | + |
| 142 | + Slope(hx, hc, h1 * u.s) |
| 143 | +end |
| 144 | + |
| 145 | +function log(u::Slope) |
| 146 | + hx = log(u.x) |
| 147 | + hc = log(u.c) |
| 148 | + |
| 149 | + i = u.x.lo - u.c.lo |
| 150 | + s = u.x.hi - u.c.hi |
| 151 | + |
| 152 | + if (i == 0 || s == 0) |
| 153 | + h1 = 1 / u.x |
| 154 | + else |
| 155 | + h1 = interval((hx.hi - hc.hi) / s, (hx.lo - hc.lo) / i) |
| 156 | + end |
| 157 | + Slope(hx, hc, h1 * u.s) |
| 158 | +end |
| 159 | + |
| 160 | +function sin(u::Slope) # Using derivative to upper bound the slope expansion for now |
| 161 | + hx = sin(u.x) |
| 162 | + hc = sin(u.c) |
| 163 | + hs = cos(u.x) |
| 164 | + Slope(hx, hc, hs) |
| 165 | +end |
| 166 | + |
| 167 | +function cos(u::Slope) # Using derivative to upper bound the slope expansion for now |
| 168 | + hx = cos(u.x) |
| 169 | + hc = cos(u.c) |
| 170 | + hs = -sin(u.x) |
| 171 | + Slope(hx, hc, hs) |
| 172 | +end |
| 173 | + |
| 174 | +function tan(u::Slope) # Using derivative to upper bound the slope expansion for now |
| 175 | + hx = tan(u.x) |
| 176 | + hc = tan(u.c) |
| 177 | + hs = (sec(u.x)) ^ 2 |
| 178 | + Slope(hx, hc, hs) |
| 179 | +end |
| 180 | + |
| 181 | +function asin(u::Slope) |
| 182 | + hx = asin(u.x) |
| 183 | + hc = asin(u.c) |
| 184 | + hs = 1 / sqrt(1 - (u.x ^ 2)) |
| 185 | + Slope(hx, hc, hs) |
| 186 | +end |
| 187 | + |
| 188 | +function acos(u::Slope) |
| 189 | + hx = acos(u.x) |
| 190 | + hc = acos(u.c) |
| 191 | + hs = -1 / sqrt(1 - (u.x ^ 2)) |
| 192 | + Slope(hx, hc, hs) |
| 193 | +end |
| 194 | + |
| 195 | +function atan(u::Slope) |
| 196 | + hx = atan(u.x) |
| 197 | + hc = atan(u.c) |
| 198 | + hs = 1 / 1 + (u.x ^ 2) |
| 199 | + Slope(hx, hc, hs) |
| 200 | +end |
0 commit comments