Object subclass: #DhbIncompleteBetaFunction instanceVariableNames: 'alpha1 alpha2 fraction inverseFraction logNorm ' classVariableNames: '' poolDictionaries: ''! Object subclass: #DhbIncompleteGammaFunction instanceVariableNames: 'alpha alphaLogGamma series fraction ' classVariableNames: '' poolDictionaries: ''! Object subclass: #DhbIterativeProcess instanceVariableNames: 'precision desiredPrecision maximumIterations result iterations ' classVariableNames: '' poolDictionaries: ''! DhbIterativeProcess subclass: #DhbFunctionalIterator instanceVariableNames: 'functionBlock relativePrecision ' classVariableNames: '' poolDictionaries: ''! DhbFunctionalIterator subclass: #DhbBisectionZeroFinder instanceVariableNames: 'positiveX negativeX ' classVariableNames: '' poolDictionaries: ''! DhbFunctionalIterator subclass: #DhbNewtonZeroFinder instanceVariableNames: 'derivativeBlock ' classVariableNames: '' poolDictionaries: ''! DhbFunctionalIterator subclass: #DhbTrapezeIntegrator instanceVariableNames: 'from to sum step ' classVariableNames: '' poolDictionaries: ''! DhbTrapezeIntegrator subclass: #DhbRombergIntegrator instanceVariableNames: 'order points interpolator ' classVariableNames: '' poolDictionaries: ''! DhbTrapezeIntegrator subclass: #DhbSimpsonIntegrator instanceVariableNames: '' classVariableNames: '' poolDictionaries: ''! DhbIterativeProcess subclass: #DhbIncompleteBetaFunctionFraction instanceVariableNames: 'x q1 q2 q3 numerator denominator alpha1 alpha2 ' classVariableNames: '' poolDictionaries: ''! DhbIterativeProcess subclass: #DhbInfiniteSeries instanceVariableNames: 'termServer ' classVariableNames: '' poolDictionaries: ''! DhbInfiniteSeries subclass: #DhbContinuedFraction instanceVariableNames: 'numerator denominator ' classVariableNames: '' poolDictionaries: ''! Object subclass: #DhbSeriesTermServer instanceVariableNames: 'x lastTerm ' classVariableNames: '' poolDictionaries: ''! DhbSeriesTermServer subclass: #DhbIncompleteBetaFractionTermServer instanceVariableNames: 'alpha1 alpha2 ' classVariableNames: '' poolDictionaries: ''! DhbSeriesTermServer subclass: #DhbIncompleteGammaFractionTermServer instanceVariableNames: 'alpha ' classVariableNames: '' poolDictionaries: ''! DhbSeriesTermServer subclass: #DhbIncompleteGammaSeriesTermServer instanceVariableNames: 'alpha sum ' classVariableNames: '' poolDictionaries: ''! SubApplication subclass: #DhbIterativeAlgorithms instanceVariableNames: '' classVariableNames: '' poolDictionaries: ''! !DhbBisectionZeroFinder publicMethods ! evaluateIteration "Perform one step of bisection. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 22/4/99 " result := ( positiveX + negativeX) * 0.5. ( functionBlock value: result) > 0 ifTrue: [ positiveX := result] ifFalse:[ negativeX := result]. ^self relativePrecision: ( positiveX - negativeX) abs! findNegativeXFrom: aNumber1 range: aNumber2 "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 22/4/99 " | n | n := 0. [ negativeX := Number random * aNumber2 + aNumber1. ( functionBlock value: negativeX) < 0 ] whileFalse: [ n := n + 0.1. n > maximumIterations ifTrue: [ self error: 'Unable to find a negative function value']. ].! findPositiveXFrom: aNumber1 range: aNumber2 "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 22/4/99 " | n | n := 0. [ positiveX := Number random * aNumber2 + aNumber1. ( functionBlock value: positiveX) > 0 ] whileFalse: [ n := n + 1. n > maximumIterations ifTrue: [ self error: 'Unable to find a positive function value']. ].! setNegativeX: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 22/4/99 " ( functionBlock value: aNumber) < 0 ifFalse:[ self error: 'Function is not negative at x = ', aNumber printString]. negativeX := aNumber.! setPositiveX: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 22/4/99 " ( functionBlock value: aNumber) > 0 ifFalse:[ self error: 'Function is not positive at x = ', aNumber printString]. positiveX := aNumber.! ! !DhbBisectionZeroFinder privateMethods ! computeInitialValues "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 22/4/99 " positiveX isNil ifTrue: [ self error: 'No positive value supplied']. negativeX isNil ifTrue: [ self error: 'No negative value supplied'].! ! !DhbContinuedFraction publicMethods ! evaluateIteration "Perform one iteration. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " | terms delta | terms := termServer termsAt: iterations. denominator := 1 / ( self limitedSmallValue: ( (terms at: 1) * denominator + (terms at: 2))). numerator := self limitedSmallValue: ( (terms at: 1) / numerator + (terms at: 2)). delta := numerator * denominator. result := result * delta. ^( delta - 1) abs! initializeIterations "Initialize the series. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " numerator := self limitedSmallValue: termServer initialTerm. denominator := 0. result := numerator! ! !DhbContinuedFraction privateMethods ! limitedSmallValue: aNumber "Private - prevent aNumber from being smaller in absolute value than a small number. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " ^aNumber abs < DhbFloatingPointMachine new smallNumber ifTrue: [ DhbFloatingPointMachine new smallNumber] ifFalse:[ aNumber]! ! !DhbFunctionalIterator class publicMethods ! function: aBlock "Convenience method to create a instance with given function block. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ^self new setFunction: aBlock; yourself! ! !DhbFunctionalIterator publicMethods ! initializeIterations "If no initial value has been defined, take 0 as the starting point (for lack of anything better). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " functionBlock isNil ifTrue: [self error: 'No function supplied']. self computeInitialValues! relativePrecision: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 27/4/99 " ^self precisionOf: aNumber relativeTo: result abs! setFunction: aBlock "Defines the function for which zeroes will be found. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ( aBlock respondsTo: #value:) ifFalse:[ self error: 'Function block must implement the method value:']. functionBlock := aBlock.! ! !DhbIncompleteBetaFractionTermServer publicMethods ! initialTerm "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " ^1! setParameter: aNumber1 second: aNumber2 "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " alpha1 := aNumber1. alpha2 := aNumber2! termsAt: anInteger "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " | n n2 | n := anInteger // 2. n2 := 2 * n. ^Array with: ( n2 < anInteger ifTrue: [x negated * (alpha1 + n) * (alpha1 + alpha2 + n) / ((alpha1 + n2) * (alpha1 + 1 + n2))] ifFalse: [x * n * (alpha2 - n) / ((alpha1 + n2) * (alpha1 - 1 + n2))]) with: 1! ! !DhbIncompleteBetaFunction class publicMethods ! shape: aNumber1 shape: aNumber2 "Create an instance of the receiver with given shape parameters. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/2/99 " ^super new initialize: aNumber1 shape: aNumber2! ! !DhbIncompleteBetaFunction publicMethods ! value: aNumber "Compute the value of the receiver for argument aNumber. Note: aNumber must be between 0 and 1 (otherwise an exception will occur) (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " | norm | aNumber = 0 ifTrue: [ ^0]. aNumber = 1 ifTrue: [ ^1]. norm := ( aNumber ln * alpha1 + ( ( 1 - aNumber) ln * alpha2) + logNorm) exp. ^( alpha1 + alpha2 + 2) * aNumber < ( alpha1 + 1) ifTrue: [ norm / ( ( self evaluateFraction: aNumber) * alpha1)] ifFalse:[ 1 - ( norm / ( ( self evaluateInverseFraction: aNumber) * alpha2))]! ! !DhbIncompleteBetaFunction privateMethods ! evaluateFraction: aNumber "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " fraction isNil ifTrue: [fraction := DhbIncompleteBetaFractionTermServer new. fraction setParameter: alpha1 second: alpha2]. fraction setArgument: aNumber. ^(DhbContinuedFraction server: fraction) desiredPrecision: DhbFloatingPointMachine new defaultNumericalPrecision; evaluate! evaluateInverseFraction: aNumber "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " inverseFraction isNil ifTrue: [inverseFraction := DhbIncompleteBetaFractionTermServer new. inverseFraction setParameter: alpha2 second: alpha1]. inverseFraction setArgument: (1 - aNumber). ^(DhbContinuedFraction server: inverseFraction) desiredPrecision: DhbFloatingPointMachine new defaultNumericalPrecision; evaluate! initialize: aNumber1 shape: aNumber2 "Private - Initialize the parameters of the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " alpha1 := aNumber1. alpha2 := aNumber2. logNorm := ( alpha1 + alpha2) logGamma - alpha1 logGamma - alpha2 logGamma. ^self! ! !DhbIncompleteBetaFunctionFraction class publicMethods ! shape: aNumber1 shape: aNumber2 "Create an instance of the receiver with given shape parameters. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/2/99 " ^super new initialize: aNumber1 shape: aNumber2! ! !DhbIncompleteBetaFunctionFraction publicMethods ! evaluateIteration "Compute and add the next term of the fraction. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " | m m2 temp | m := iterations + 1. m2 := m * 2. temp := m * (alpha2 - m) * x / ((q3 + m2) * (alpha1 + m2)). denominator := self limitedSmallValue: ( denominator * temp + 1). numerator := self limitedSmallValue: ( temp / numerator + 1). denominator := 1 / denominator. result := result * numerator * denominator. temp := (alpha1 + m) negated * (q1 + m) * x / ((q2 + m2) * (alpha1 + m2)). denominator := self limitedSmallValue: ( denominator * temp + 1). numerator := self limitedSmallValue: ( temp / numerator + 1). denominator := 1 / denominator. temp := numerator * denominator. result := result * temp. ^(temp - 1) abs! initializeIterations "Initialize the iterations (subclasses must write their own method and call this one last). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " numerator := 1. denominator := 1 / (self limitedSmallValue: 1 - (q1 * x / q2)). result := denominator! setArgument: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " x := aNumber.! ! !DhbIncompleteBetaFunctionFraction privateMethods ! initialize: aNumber1 shape: aNumber2 "Private - Initialize the parameters of the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " alpha1 := aNumber1. alpha2 := aNumber2. q1 := alpha1 + alpha2. q2 := alpha1 + 1. q3 := alpha1 - 1. ^self! ! !DhbIncompleteGammaFractionTermServer publicMethods ! initialTerm "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " lastTerm := x - alpha + 1. ^lastTerm! setParameter: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " alpha := aNumber asFloat! termsAt: anInteger "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " lastTerm := lastTerm + 2. ^Array with: (alpha - anInteger) * anInteger with: lastTerm! ! !DhbIncompleteGammaFunction class publicMethods ! shape: aNumber "Defines a new instance of the receiver with paramater aNumber (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " ^super new initialize: aNumber! ! !DhbIncompleteGammaFunction publicMethods ! value: aNumber "Compute the value of the receiver for argument aNumber. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " | x norm | aNumber = 0 ifTrue: [ ^0]. x := aNumber asFloat. norm := [ ( x ln * alpha - x - alphaLogGamma) exp] when: ExAll do: [ :signal | signal exitWith: nil]. norm isNil ifTrue: [ ^1]. ^x - 1 < alpha ifTrue: [ ( self evaluateSeries: x) * norm] ifFalse:[ 1 - ( norm / ( self evaluateFraction: x))]! ! !DhbIncompleteGammaFunction privateMethods ! evaluateFraction: aNumber "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " fraction isNil ifTrue: [fraction := DhbIncompleteGammaFractionTermServer new. fraction setParameter: alpha]. fraction setArgument: aNumber. ^(DhbContinuedFraction server: fraction) desiredPrecision: DhbFloatingPointMachine new defaultNumericalPrecision; evaluate! evaluateSeries: aNumber "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " series isNil ifTrue: [ series := DhbIncompleteGammaSeriesTermServer new. series setParameter: alpha. ]. series setArgument: aNumber. ^(DhbInfiniteSeries server: series) desiredPrecision: DhbFloatingPointMachine new defaultNumericalPrecision; evaluate! initialize: aNumber "Private - Defines the parameter alpha of the receiver (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 10/3/99 " alpha := aNumber asFloat. alphaLogGamma := alpha logGamma. ^self! ! !DhbIncompleteGammaSeriesTermServer publicMethods ! initialTerm "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " lastTerm := 1 / alpha. sum := alpha. ^lastTerm! setParameter: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " alpha := aNumber asFloat! termAt: anInteger "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " sum := sum + 1. lastTerm := lastTerm * x / sum. ^lastTerm! ! !DhbInfiniteSeries class publicMethods ! server: aTermServer "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " ^self new initialize: aTermServer! ! !DhbInfiniteSeries publicMethods ! evaluateIteration "Perform one iteration. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " | delta | delta := termServer termAt: iterations. result := result + delta. ^self precisionOf: delta abs relativeTo: result abs! initializeIterations "Initialize the series. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " result := termServer initialTerm! ! !DhbInfiniteSeries privateMethods ! initialize: aTermServer "Private - Assigns the object responsible to compute each term. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " termServer := aTermServer. ^self! ! !DhbIterativeProcess class publicMethods ! new "Create an instance of the class. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^super new initialize! ! !DhbIterativeProcess class privateMethods ! defaultMaximumIterations "Private - Answers the default maximum number of iterations for newly created instances. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^50! defaultPrecision "Private - Answers the default precision for newly created instances. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^DhbFloatingPointMachine new defaultNumericalPrecision! ! !DhbIterativeProcess publicMethods ! desiredPrecision: aNumber "Defines the desired precision for the result. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " aNumber > 0 ifFalse: [ ^self error: 'Illegal precision: ', aNumber printString]. desiredPrecision := aNumber.! evaluate "Perform the iteration until either the desired precision is attained or the number of iterations exceeds the maximum. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " iterations := 0. self initializeIterations. [iterations := iterations + 1. precision := self evaluateIteration. self hasConverged or: [iterations >= maximumIterations]] whileFalse: []. self finalizeIterations. ^self result! evaluateIteration "Dummy method (must be implemented by subclass). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^self subclassResponsibility! finalizeIterations "Perform cleanup operation if needed (must be implemented by subclass). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 9/3/99 " ! hasConverged " (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 20/4/99 " ^precision <= desiredPrecision! initializeIterations "Initialize the iterations (must be implemented by subclass when needed). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ! iterations "Answers the number of iterations performed. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^iterations! maximumIterations: anInteger "Defines the maximum number of iterations. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ( anInteger isInteger and: [ anInteger > 1]) ifFalse: [ ^self error: 'Invalid maximum number of iteration: ', anInteger printString]. maximumIterations := anInteger.! precision "Answer the attained precision for the result. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ^precision! precisionOf: aNumber1 relativeTo: aNumber2 "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 3/5/99 " ^aNumber2 > DhbFloatingPointMachine new defaultNumericalPrecision ifTrue: [ aNumber1 / aNumber2] ifFalse:[ aNumber1]! result "Answer the result of the iterations (if any) (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^result! ! !DhbIterativeProcess privateMethods ! initialize "Private - initialize the parameters of the receiver with default values. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " desiredPrecision := self class defaultPrecision. maximumIterations := self class defaultMaximumIterations. ^self! ! !DhbNewtonZeroFinder class publicMethods ! function: aBlock1 derivative: aBlock2 "Convenience method to create a instance with given function block. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ^(self new) setFunction: aBlock1; setDerivative: aBlock2; yourself! ! !DhbNewtonZeroFinder publicMethods ! evaluateIteration "Compute one step of Newton's zero finding method. Answers the estimated precision. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " | delta | delta := ( functionBlock value: result) / ( derivativeBlock value: result). result := result - delta. ^self relativePrecision: delta abs! initialValue: aNumber "Define the initial value for the iterations. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " result := aNumber.! setDerivative: aBlock "Defines the derivative of the function for which zeroes will be found. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " | x | ( aBlock respondsTo: #value:) ifFalse:[ self error: 'Derivative block must implement the method value:']. x := result ifNil: [ Number random] ifNot: [ :base | base + Number random]. ( ( aBlock value: x) relativelyEqualsTo: (self defaultDerivativeBlock value: x) upTo: 0.0001) ifFalse:[ self error: 'Supplied derivative is not correct']. derivativeBlock := aBlock.! setFunction: aBlock "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " super setFunction: aBlock. derivativeBlock := nil.! ! !DhbNewtonZeroFinder privateMethods ! computeInitialValues "Private - If no derivative has been defined, take an ad-hoc definition. If no initial value has been defined, take 0 as the starting point (for lack of anything better). (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " | n | result isNil ifTrue: [ result := 0]. derivativeBlock isNil ifTrue: [ derivativeBlock := self defaultDerivativeBlock]. n := 0. [ (derivativeBlock value: result) equalsTo: 0] whileTrue: [ n := n + 1. n > maximumIterations ifTrue: [ self error: 'Function''s derivative seems to be zero everywhere']. result := Number random + result].! defaultDerivativeBlock "Private - Answers a block computing the function's derivative by approximation. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ^[ :x | 5000 * ( ( functionBlock value: (x + 0.0001)) - ( functionBlock value: (x - 0.0001)))]! ! !DhbPolynomial publicMethods ! deflatedAt: aNumber "Answers a new polynomial quotient of the receiver with polynomial (X-aNumber) (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 20/4/99 " | remainder next newCoefficients| remainder := 0. newCoefficients := coefficients collect: [ :each | next := remainder. remainder := remainder * aNumber + each. next]. ^self class new: ( newCoefficients copyFrom: 2 to: newCoefficients size) reverse! roots "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 20/4/99 " ^self roots: DhbFloatingPointMachine new defaultNumericalPrecision ! roots: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 20/4/99 " | pol roots x rootFinder | rootFinder := DhbNewtonZeroFinder new. rootFinder desiredPrecision: aNumber. pol := self class new: ( coefficients reverse collect: [ :each | each asFloat]). roots := OrderedCollection new: self degree. [ rootFinder setFunction: pol; setDerivative: pol derivative. x := rootFinder evaluate. rootFinder hasConverged ] whileTrue: [ roots add: x. pol := pol deflatedAt: x. pol degree > 0 ifFalse: [ ^roots]. ]. ^roots ! ! !DhbRombergIntegrator class privateMethods ! defaultOrder "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " ^5! ! !DhbRombergIntegrator publicMethods ! evaluateIteration "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " | interpolation | points addLast: (points last x * 0.25) @ self higherOrderSum. points size < order ifTrue: [ ^1]. interpolation := interpolator valueAndError: 0. points removeFirst. result := interpolation at: 1. ^self relativePrecision: ( interpolation at: 2) abs! order: anInteger "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " anInteger < 2 ifTrue: [ self error: 'Order for Romberg integration must be larger than 1']. order := anInteger.! ! !DhbRombergIntegrator privateMethods ! computeInitialValues "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " super computeInitialValues. points := OrderedCollection new: order. interpolator := DhbNevilleInterpolator points: points. points add: 1 @ sum. ! initialize "Private - initialize the parameters of the receiver with default values. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " order := self class defaultOrder. ^super initialize! ! !DhbSeriesTermServer publicMethods ! setArgument: aNumber "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 15/3/99 " x := aNumber asFloat.! ! !DhbSimpsonIntegrator publicMethods ! evaluateIteration "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 27/4/99 " | oldResult oldSum | iterations < 2 ifTrue: [ self higherOrderSum. ^1 ]. oldResult := result. oldSum := sum. result := (self higherOrderSum * 4 - oldSum) / 3. ^self relativePrecision: ( result - oldResult) abs! ! !DhbTrapezeIntegrator class publicMethods ! new: aBlock from: aNumber1 to: aNumber2 "Create an new instance with given parameters. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ^super new initialize: aBlock from: aNumber1 to: aNumber2! ! !DhbTrapezeIntegrator class privateMethods ! defaultMaximumIterations "Private - Answers the default maximum number of iterations for newly created instances. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 6/1/99 " ^13! new "Private - Block the constructor method for this class. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " ^self error: 'Method new:from:to: must be used'! ! !DhbTrapezeIntegrator publicMethods ! evaluateIteration "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " | oldResult | oldResult := result. result := self higherOrderSum. ^self relativePrecision: ( result - oldResult) abs! from: aNumber1 to: aNumber2 "(c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " from := aNumber1. to := aNumber2.! ! !DhbTrapezeIntegrator privateMethods ! computeInitialValues "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " step := to - from. sum := ( ( functionBlock value: from) + ( functionBlock value: to)) * step /2. result := sum.! higherOrderSum "Private - (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 26/4/99 " | x newSum | x := step / 2 + from. newSum := 0. [ x < to ] whileTrue: [ newSum := ( functionBlock value: x) + newSum. x := x + step. ]. sum := ( step * newSum + sum) / 2. step := step / 2. ^sum ! initialize: aBlock from: aNumber1 to: aNumber2 "Private - Initialize the parameters of the receiver. (c) Copyrights Didier BESSET, 1999, all rights reserved. Initial code: 7/1/99 " functionBlock := aBlock. self from: aNumber1 to: aNumber2. ^self ! ! DhbIncompleteBetaFunction initializeAfterLoad! DhbIncompleteGammaFunction initializeAfterLoad! DhbIterativeProcess initializeAfterLoad! DhbFunctionalIterator initializeAfterLoad! DhbBisectionZeroFinder initializeAfterLoad! DhbNewtonZeroFinder initializeAfterLoad! DhbTrapezeIntegrator initializeAfterLoad! DhbRombergIntegrator initializeAfterLoad! DhbSimpsonIntegrator initializeAfterLoad! DhbIncompleteBetaFunctionFraction initializeAfterLoad! DhbInfiniteSeries initializeAfterLoad! DhbContinuedFraction initializeAfterLoad! DhbSeriesTermServer initializeAfterLoad! DhbIncompleteBetaFractionTermServer initializeAfterLoad! DhbIncompleteGammaFractionTermServer initializeAfterLoad! DhbIncompleteGammaSeriesTermServer initializeAfterLoad! DhbIterativeAlgorithms initializeAfterLoad!