Skip to content

Commit f8bffee

Browse files
kylecorry31kylecorry31
committed
Add numeric differentiation and polynomials
1 parent 5cb79eb commit f8bffee

File tree

6 files changed

+484
-62
lines changed

6 files changed

+484
-62
lines changed

‎src/main/kotlin/com/kylecorry/sol/math/algebra/Polynomial.kt‎

Lines changed: 143 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,143 @@
1+
package com.kylecorry.sol.math.algebra
2+
3+
import com.kylecorry.sol.math.SolMath
4+
import com.kylecorry.sol.math.sumOfFloat
5+
6+
data class PolynomialTerm(
7+
val coefficient: Float,
8+
val exponent: Int
9+
)
10+
11+
class Polynomial(terms: List<PolynomialTerm>) {
12+
13+
val terms = terms
14+
.groupBy { it.exponent }
15+
.map { (exponent, groupedTerms) ->
16+
val coefficient = groupedTerms.sumOfFloat { it.coefficient }
17+
PolynomialTerm(coefficient, exponent)
18+
}
19+
.filter { !SolMath.isZero(it.coefficient) }
20+
.sortedBy { -it.exponent }
21+
22+
fun evaluate(x: Float): Float {
23+
return terms.sumOfFloat { it.coefficient * SolMath.power(x, it.exponent) }
24+
}
25+
26+
fun derivative(): Polynomial {
27+
val derivedTerms = terms.mapNotNull {
28+
if (it.exponent == 0) {
29+
null
30+
} else {
31+
PolynomialTerm(it.coefficient * it.exponent, it.exponent - 1)
32+
}
33+
}
34+
return Polynomial(derivedTerms)
35+
}
36+
37+
fun integral(c: Float = 0f): Polynomial {
38+
val integratedTerms = terms.map {
39+
PolynomialTerm(
40+
coefficient = it.coefficient / (it.exponent + 1),
41+
exponent = it.exponent + 1
42+
)
43+
}
44+
return Polynomial(integratedTerms + listOf(PolynomialTerm(c, 0)))
45+
}
46+
47+
override fun equals(other: Any?): Boolean {
48+
if (this === other) return true
49+
if (other !is Polynomial) return false
50+
51+
if (terms.size != other.terms.size) return false
52+
for (i in terms.indices) {
53+
if (terms[i] != other.terms[i]) return false
54+
}
55+
56+
return true
57+
}
58+
59+
override fun hashCode(): Int {
60+
return terms.hashCode()
61+
}
62+
63+
override fun toString(): String {
64+
if (terms.isEmpty()) {
65+
return "0"
66+
}
67+
return terms.joinToString(" + ") { term ->
68+
val coefficientString = when {
69+
term.exponent == 0 -> "${term.coefficient}"
70+
term.coefficient == 1f -> ""
71+
term.coefficient == -1f -> "-"
72+
else -> "${term.coefficient}"
73+
}.removeSuffix(".0")
74+
val exponentStr = when (term.exponent) {
75+
0 -> ""
76+
1 -> "x"
77+
else -> "x^${term.exponent}"
78+
}
79+
"$coefficientString$exponentStr"
80+
}.replace("+ -", "- ")
81+
}
82+
83+
companion object {
84+
fun of(vararg terms: PolynomialTerm): Polynomial {
85+
return Polynomial(terms.toList())
86+
}
87+
88+
/**
89+
* Creates a polynomial from coefficients.
90+
* The index of the coefficient corresponds to the exponent.
91+
*/
92+
fun fromCoefficients(vararg coefficients: Float): Polynomial {
93+
val terms = coefficients.mapIndexed { index, coefficient ->
94+
PolynomialTerm(coefficient, index)
95+
}
96+
return Polynomial(terms)
97+
}
98+
99+
/**
100+
* Creates a polynomial from a string equation.
101+
* Example inputs: "x^2", "x + 4x^2 + 5", "-2x", "x - 2x^3"
102+
*/
103+
fun of(equation: String): Polynomial {
104+
val terms = equation
105+
.replace(" ", "")
106+
.replace("-", "+-")
107+
.replace("^", "")
108+
.split("+")
109+
.filter { it.isNotEmpty() }
110+
.map { termString ->
111+
if ("x" !in termString) {
112+
return@map PolynomialTerm(termString.toFloat(), 0)
113+
}
114+
val parts = termString.split("x")
115+
val coefficientPart = parts[0]
116+
.replace("(", "")
117+
.replace(")", "")
118+
val coefficient = when {
119+
parts[0] == "" -> 1f
120+
parts[0] == "-" -> -1f
121+
parts[0].contains("/") -> {
122+
val fractionParts = coefficientPart.split("/")
123+
if (fractionParts.size != 2) {
124+
throw IllegalArgumentException("Invalid fraction in term: $termString")
125+
}
126+
val numerator = fractionParts[0].toFloat()
127+
val denominator = fractionParts[1].toFloat()
128+
numerator / denominator
129+
}
130+
131+
else -> parts[0].toFloat()
132+
}
133+
val exponent = if (parts.size > 1 && parts[1].isNotEmpty()) {
134+
parts[1].toInt()
135+
} else {
136+
1
137+
}
138+
PolynomialTerm(coefficient, exponent)
139+
}
140+
return Polynomial(terms)
141+
}
142+
}
143+
}

‎src/main/kotlin/com/kylecorry/sol/math/calculus/Calculus.kt‎

Lines changed: 17 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
11
package com.kylecorry.sol.math.calculus
22

3-
import com.kylecorry.sol.math.SolMath
43
import com.kylecorry.sol.math.Vector2
54
import com.kylecorry.sol.math.algebra.LinearEquation
5+
import com.kylecorry.sol.math.algebra.Polynomial
66
import com.kylecorry.sol.math.algebra.QuadraticEquation
77
import com.kylecorry.sol.shared.Guards
88
import kotlin.math.abs
@@ -11,6 +11,10 @@ import kotlin.math.min
1111

1212
object Calculus {
1313

14+
fun derivative(polynomial: Polynomial): Polynomial {
15+
return polynomial.derivative()
16+
}
17+
1418
fun derivative(equation: QuadraticEquation): LinearEquation {
1519
return LinearEquation(equation.a * 2, equation.b)
1620
}
@@ -22,23 +26,15 @@ object Calculus {
2226
/**
2327
* Get the derivative of samples.
2428
* This assumes values is sorted by x (increasing).
25-
* The resulting derivative will have one value less (first x value dropped)
29+
* @param values The samples to differentiate
30+
* @param finiteDifferenceOrder The order of the finite difference to use (1, 2, or 3). Default is 1.
31+
* @return The derivative samples
2632
*/
27-
fun derivative(values: List<Vector2>): List<Vector2> {
28-
val derivative = mutableListOf<Vector2>()
29-
for (j in 0 until values.size - 1) {
30-
val x = values[j + 1].x
31-
val dx = values[j + 1].x - values[j].x
32-
if (SolMath.isZero(dx)) {
33-
derivative.add(Vector2(x, 0f))
34-
}
35-
36-
val dy = (values[j + 1].y - values[j].y) / dx
37-
derivative.add(Vector2(x, dy))
38-
}
39-
return derivative
33+
fun derivative(values: List<Vector2>, finiteDifferenceOrder: Int = 1): List<Vector2> {
34+
return NumericDifferentiation.finiteDifference(values, finiteDifferenceOrder)
4035
}
4136

37+
// TODO: Use central difference
4238
fun derivative(
4339
x: Double,
4440
step: Double = 0.0001,
@@ -48,6 +44,7 @@ object Calculus {
4844
return (fn(x + step) - current) / step
4945
}
5046

47+
// TODO: Use central difference
5148
fun derivative(
5249
x: Double,
5350
y: Double,
@@ -60,6 +57,7 @@ object Calculus {
6057
return xGrad to yGrad
6158
}
6259

60+
// TODO: RK4 solver?
6361
fun integral(
6462
startX: Double,
6563
endX: Double,
@@ -101,6 +99,10 @@ object Calculus {
10199
return multiplier * total
102100
}
103101

102+
fun integral(polynomial: Polynomial, c: Float = 0f): Polynomial {
103+
return polynomial.integral(c)
104+
}
105+
104106
/**
105107
* Calculate the root of the provided function using Newton's Method
106108
*/
@@ -120,5 +122,4 @@ object Calculus {
120122
return x
121123
}
122124

123-
124125
}

‎src/main/kotlin/com/kylecorry/sol/math/calculus/NumericDifferentiation.kt‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
1+
package com.kylecorry.sol.math.calculus
2+
3+
import com.kylecorry.sol.math.SolMath
4+
import com.kylecorry.sol.math.Vector2
5+
import kotlin.math.max
6+
import kotlin.math.min
7+
8+
internal object NumericDifferentiation {
9+
10+
private val FORWARD_STENCILS = mapOf(
11+
1 to Stencil(
12+
intArrayOf(0, 1),
13+
floatArrayOf(-1f, 1f)
14+
),
15+
2 to Stencil(
16+
intArrayOf(0, 1, 2),
17+
floatArrayOf(-3f / 2f, 2f, -1f / 2f)
18+
),
19+
3 to Stencil(
20+
intArrayOf(0, 1, 2, 3),
21+
floatArrayOf(-11f / 6f, 3f, -3f / 2f, 1f / 3f)
22+
)
23+
)
24+
25+
private val CENTRAL_STENCILS = mapOf(
26+
1 to Stencil(
27+
intArrayOf(-1, 1),
28+
floatArrayOf(-0.5f, 0.5f)
29+
),
30+
2 to Stencil(
31+
intArrayOf(-2, -1, 1, 2),
32+
floatArrayOf(1f / 12f, -2f / 3f, 2f / 3f, -1f / 12f)
33+
),
34+
3 to Stencil(
35+
intArrayOf(-3, -2, -1, 1, 2, 3),
36+
floatArrayOf(-1f / 60f, 3f / 20f, -3f / 4f, 3f / 4f, -3f / 20f, 1f / 60f)
37+
)
38+
)
39+
40+
private val BACKWARD_STENCILS = mapOf(
41+
1 to Stencil(
42+
intArrayOf(-1, 0),
43+
floatArrayOf(-1f, 1f)
44+
),
45+
2 to Stencil(
46+
intArrayOf(-2, -1, 0),
47+
floatArrayOf(1f / 2f, -2f, 3f / 2f)
48+
),
49+
3 to Stencil(
50+
intArrayOf(-3, -2, -1, 0),
51+
floatArrayOf(-1f / 3f, 3f / 2f, -3f, 11f / 6f)
52+
)
53+
)
54+
55+
private fun stencilDifference(
56+
values: List<Vector2>,
57+
i: Int,
58+
stencil: Stencil
59+
): Vector2 {
60+
// TODO: dx multiplier
61+
var derivative = 0f
62+
var minX = Float.MAX_VALUE
63+
var maxX = -Float.MAX_VALUE
64+
val x = values[i].x
65+
val range = stencil.indices.last() - stencil.indices.first()
66+
for (j in stencil.indices.indices) {
67+
val index = i + stencil.indices[j]
68+
if (index < 0 || index >= values.size) {
69+
continue
70+
}
71+
derivative += values[index].y * stencil.coefficients[j]
72+
minX = min(minX, values[index].x)
73+
maxX = max(maxX, values[index].x)
74+
}
75+
val dx = maxX - minX
76+
if (SolMath.isZero(dx)) {
77+
return Vector2(x, 0f)
78+
}
79+
// Coefficients are based on the multiple of dx (space between points, which is why it is count - 1)
80+
derivative *= range
81+
return Vector2(x, derivative / dx)
82+
}
83+
84+
fun finiteDifference(values: List<Vector2>, order: Int = 1): List<Vector2> {
85+
if (values.size < 2) {
86+
return emptyList()
87+
}
88+
val actualOrder = order.coerceIn(1, 3)
89+
90+
val derivative = mutableListOf<Vector2>()
91+
// First point needs forward difference
92+
derivative.add(stencilDifference(values, 0, FORWARD_STENCILS[actualOrder]!!))
93+
for (j in 1 until values.size - 1) {
94+
val neighbors = min(j, values.size - 1 - j)
95+
val centerOrder = min(actualOrder, neighbors)
96+
derivative.add(stencilDifference(values, j, CENTRAL_STENCILS[centerOrder]!!))
97+
}
98+
// Last point needs backward difference
99+
derivative.add(stencilDifference(values, values.lastIndex, BACKWARD_STENCILS[actualOrder]!!))
100+
return derivative
101+
}
102+
103+
// Second derivative
104+
// Central: (f(t + dt) - 2f(t) + f(t - dt)) / dt^2
105+
106+
class Stencil(
107+
val indices: IntArray,
108+
val coefficients: FloatArray
109+
)
110+
111+
}

0 commit comments

Comments
 (0)