Use split interval for integration.

This commit is contained in:
Alexander Nozik 2021-04-18 19:43:03 +03:00
parent 6c215abf13
commit 07e39a068d
5 changed files with 56 additions and 18 deletions

View File

@ -22,8 +22,8 @@ fun main() {
return DoubleBuffer(x.size) { i -> return DoubleBuffer(x.size) { i ->
val h = sigma[i] / 5 val h = sigma[i] / 5
val dVector = DoubleBuffer(x.size) { if (it == i) h else 0.0 } val dVector = DoubleBuffer(x.size) { if (it == i) h else 0.0 }
val f1 = invoke(x + dVector / 2) val f1 = this(x + dVector / 2)
val f0 = invoke(x - dVector / 2) val f0 = this(x - dVector / 2)
(f1 - f0) / h (f1 - f0) / h
} }
} }

View File

@ -28,7 +28,7 @@ public class CMIntegrator(
val integrator = integratorBuilder(integrand) val integrator = integratorBuilder(integrand)
val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls val maxCalls = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: defaultMaxCalls
val remainingCalls = maxCalls - integrand.calls val remainingCalls = maxCalls - integrand.calls
val range = integrand.getFeature<IntegrationRange<Double>>()?.range val range = integrand.getFeature<IntegrationRange>()?.range
?: error("Integration range is not provided") ?: error("Integration range is not provided")
val res = integrator.integrate(remainingCalls, integrand.function, range.start, range.endInclusive) val res = integrator.integrate(remainingCalls, integrand.function, range.start, range.endInclusive)

View File

@ -17,7 +17,7 @@ public class GaussRuleIntegrator(
) : UnivariateIntegrator<Double> { ) : UnivariateIntegrator<Double> {
override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> { override fun integrate(integrand: UnivariateIntegrand<Double>): UnivariateIntegrand<Double> {
val range = integrand.getFeature<IntegrationRange<Double>>()?.range val range = integrand.getFeature<IntegrationRange>()?.range
?: error("Integration range is not provided") ?: error("Integration range is not provided")
val integrator: GaussIntegrator = getIntegrator(range) val integrator: GaussIntegrator = getIntegrator(range)
//TODO check performance //TODO check performance

View File

@ -8,6 +8,13 @@ import space.kscience.kmath.misc.UnstableKMathAPI
import space.kscience.kmath.operations.Field import space.kscience.kmath.operations.Field
import space.kscience.kmath.structures.* import space.kscience.kmath.structures.*
/**
* Set of univariate integration ranges. First components correspond to ranges themselves, second components to number of
* integration nodes per range
*/
public class UnivariateIntegrandRanges(public val ranges: List<Pair<ClosedRange<Double>, Int>>) : IntegrandFeature {
public constructor(vararg pairs: Pair<ClosedRange<Double>, Int>) : this(pairs.toList())
}
/** /**
* A simple one-pass integrator based on Gauss rule * A simple one-pass integrator based on Gauss rule
@ -18,9 +25,29 @@ public class GaussIntegrator<T : Any>(
private fun buildRule(integrand: UnivariateIntegrand<T>): Pair<Buffer<Double>, Buffer<Double>> { private fun buildRule(integrand: UnivariateIntegrand<T>): Pair<Buffer<Double>, Buffer<Double>> {
val factory = integrand.getFeature<GaussIntegratorRuleFactory>() ?: GaussLegendreRuleFactory val factory = integrand.getFeature<GaussIntegratorRuleFactory>() ?: GaussLegendreRuleFactory
val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100 val predefinedRanges = integrand.getFeature<UnivariateIntegrandRanges>()
val range = integrand.getFeature<IntegrationRange<Double>>()?.range ?: 0.0..1.0 if (predefinedRanges == null || predefinedRanges.ranges.isEmpty()) {
return factory.build(numPoints, range) val numPoints = integrand.getFeature<IntegrandMaxCalls>()?.maxCalls ?: 100
val range = integrand.getFeature<IntegrationRange>()?.range ?: 0.0..1.0
return factory.build(numPoints, range)
} else {
val ranges = predefinedRanges.ranges
var counter = 0
val length = ranges.sumOf { it.second }
val pointsArray = DoubleArray(length)
val weightsArray = DoubleArray(length)
for (range in ranges) {
val rule = factory.build(range.second, range.first)
repeat(rule.first.size) { i ->
pointsArray[counter] = rule.first[i]
weightsArray[counter] = rule.second[i]
counter++
}
}
return pointsArray.asBuffer() to weightsArray.asBuffer()
}
} }
override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) { override fun integrate(integrand: UnivariateIntegrand<T>): UnivariateIntegrand<T> = with(algebra) {
@ -49,7 +76,8 @@ public class GaussIntegrator<T : Any>(
* Following features are evaluated: * Following features are evaluated:
* * [GaussIntegratorRuleFactory] - A factory for computing the Gauss integration rule. By default uses [GaussLegendreRuleFactory] * * [GaussIntegratorRuleFactory] - A factory for computing the Gauss integration rule. By default uses [GaussLegendreRuleFactory]
* * [IntegrationRange] - the univariate range of integration. By default uses 0..1 interval. * * [IntegrationRange] - the univariate range of integration. By default uses 0..1 interval.
* * [IntegrandMaxCalls] - the maximum number of function calls during integration. For non-iterative rules, always uses the maximum number of points. By default uses 100 points. * * [IntegrandMaxCalls] - the maximum number of function calls during integration. For non-iterative rules, always uses the maximum number of points. By default uses 10 points.
* * [UnivariateIntegrandRanges] - Set of ranges and number of points per range. Defaults to given [IntegrationRange] and [IntegrandMaxCalls]
*/ */
@UnstableKMathAPI @UnstableKMathAPI
public fun <T : Any> Field<T>.integrate( public fun <T : Any> Field<T>.integrate(
@ -64,15 +92,25 @@ public fun <T : Any> Field<T>.integrate(
@UnstableKMathAPI @UnstableKMathAPI
public fun <T : Any> Field<T>.integrate( public fun <T : Any> Field<T>.integrate(
range: ClosedRange<Double>, range: ClosedRange<Double>,
numPoints: Int = 100, order: Int = 10,
intervals: Int = 10,
vararg features: IntegrandFeature, vararg features: IntegrandFeature,
function: (Double) -> T, function: (Double) -> T,
): UnivariateIntegrand<T> = GaussIntegrator(this).integrate( ): UnivariateIntegrand<T> {
UnivariateIntegrand( require(range.endInclusive > range.start) { "The range upper bound should be higher than lower bound" }
function, require(order > 1) { "The order of polynomial must be more than 1" }
IntegrationRange(range), require(intervals > 0) { "Number of intervals must be positive" }
GaussLegendreRuleFactory, val rangeSize = (range.endInclusive - range.start) / intervals
IntegrandMaxCalls(numPoints), val ranges = UnivariateIntegrandRanges(
*features (0 until intervals).map { i -> (rangeSize * i)..(rangeSize * i + 1) to order }
) )
) return GaussIntegrator(this).integrate(
UnivariateIntegrand(
function,
IntegrationRange(range),
GaussLegendreRuleFactory,
ranges,
*features
)
)
}

View File

@ -33,7 +33,7 @@ public fun <T : Any> UnivariateIntegrand(
public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>> public typealias UnivariateIntegrator<T> = Integrator<UnivariateIntegrand<T>>
@JvmInline @JvmInline
public value class IntegrationRange<T : Comparable<T>>(public val range: ClosedRange<T>) : IntegrandFeature public value class IntegrationRange(public val range: ClosedRange<Double>) : IntegrandFeature
public val <T : Any> UnivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value public val <T : Any> UnivariateIntegrand<T>.value: T? get() = getFeature<IntegrandValue<T>>()?.value