dichotomie.ts 16.59 KiB
import { BoolIdentity, Debug } from "./internal_modules";
import { Nub } from "./internal_modules";
import { ParamDefinition } from "./internal_modules";
import { ParamDomain, ParamDomainValue } from "./internal_modules";
import { SessionSettings } from "./internal_modules";
import { Interval } from "./internal_modules";
import { Message, MessageCode } from "./internal_modules";
import { Result } from "./internal_modules";
import { SearchInterval } from "./internal_modules";
/**
 * calcul par dichotomie
 * principe : y=f(x1,x2,...), on connait une valeur de y,
 * on cherche par ex le x2 correspondant, les autres xi étant fixés.
 * la méthode Equation() calcule analytiquement y=f(xi)
export class Dichotomie extends Debug {
    /**
     * définition de la variable de la fonction
    private _paramX: ParamDefinition;
    /**
     * nombre d'étapes de recherche de l'intervalle de départ
    private _startIntervalMaxSteps = 50;
    /**
     * nom du paramètre calculé analytiquement par Equation()
    private analyticalSymbol: string;
    /**
     * Targetted value calculated by the function
    private _target: number;
    private _func: () => number;
    /**
     * si > 0, un calcul dichotomique est en cours
    private static _inDicho: number = 0;
    /**
     * Construction de la classe.
     * @param nub Noeud de calcul contenant la méthode de calcul Equation
     * @param sVarCalc Nom de la variable à calculer
    constructor(private nub: Nub, private sVarCalc: string, dbg: boolean = false, func?: () => number) {
        super(dbg);
        this._paramX = this.nub.getParameter(this.sVarCalc);
        if (func !== undefined) {
            this._func = func;
        } else {
            this._func = this.nub.EquationFirstAnalyticalParameter;
            this.analyticalSymbol = this.nub.firstAnalyticalPrmSymbol;
    private set vX(vCalc: number) {
        this._paramX.v = vCalc;
    /**
     * étapes de recherche de l'intervalle de départ
    get startIntervalMaxSteps() {
        return this._startIntervalMaxSteps;
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
set startIntervalMaxSteps(n: number) { this._startIntervalMaxSteps = n; } public static get inDico(): boolean { return Dichotomie._inDicho > 0; } public CalculX(x: number): number { this.vX = x; return this.Calcul(); } /** * trouve la valeur x pour y=f(x) avec un y donné * @param rTarget valeur de y * @param rTol tolérance pour l'arrêt de la recherche * @param rInit valeur initiale approximative de x */ public Dichotomie(rTarget: number, rTol: number, rInit: number): Result { Dichotomie._inDicho++; this._target = rTarget; // console.log("-----"); // for (let x = 0; x <= 1; x += 0.1) // console.log(this.CalculX(x).vCalc); // console.log("-----"); // recherche de l'intervalle de départ this.debug( "Valeur cible de " + this.analyticalSymbol + "=" + rTarget + ", valeur initiale de " + this.sVarCalc + "=" + rInit ); try { const r = this.getStartInterval(rTarget, rInit); if (r.ok) { const inter: SearchInterval = r.intSearch; // Recherche du zero par la méthode de Brent const s = this.uniroot( this.CalcZero, this, inter.min, inter.max, rTol ); if (s.found) { this.debug( `${this.analyticalSymbol}(${this.sVarCalc}=${s.value}) = ${rTarget}` ); Dichotomie._inDicho--; return new Result(s.value); } else { this.debug( "Non convergence" ); Dichotomie._inDicho--; return new Result(new Message(MessageCode.ERROR_DICHO_CONVERGE, { lastApproximation: s.value }), this.nub); } } else { Dichotomie._inDicho--; return new Result(r.res); } } catch (e) { // un appel à Calcul() a généré une erreur Dichotomie._inDicho--;
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
return this.nub.result; } } public debug(s: string) { if (this.DBG) { super.debug("Dichotomie: " + s); } } /** * Calcul de l'équation analytique. * @note Wrapper vers this.nub.Equation pour simplifier le code. * On utilise la première variable du tableau des variables pouvant être calculée analytiquement * Il faudra s'assurer que cette première variable correspond à la méthode de calcul la plus rapide */ private Calcul(): number { // return this.nub.Equation(this.analyticalSymbol).vCalc; return this._func.call(this.nub); } /** * Détermine si une fonction est croissante ou décroissante. * @param x point auquel on calcul la dérivée * @param dom domaine de définition de la variable */ private isIncreasingFunction(x: number, dom: Interval): boolean { let epsilon = 1e-8; for (let i = 0; i < 20; i++) { const bounds = new Interval(x - epsilon, x + epsilon); bounds.setInterval(bounds.intersect(dom)); // au cas où l'on sorte du domaine de la variable de la fonction const y1 = this.CalculX(bounds.min); const y2 = this.CalculX(bounds.max); if (Math.abs(y2 - y1) > 1E-6) return y2 > y1; epsilon *= 10; } return true; } /** * recherche l'intervalle contenant la valeur cible * @param rTarget valeur cible * @param intSearch intervalle de départ * @param intMax intervalle maximum (domaine de définition de la variable) */ private searchTarget(rTarget: number, intSearch: SearchInterval, intMax: Interval) { this.debug(`searchTarget debut: Target ${this.analyticalSymbol}=${rTarget} dans` + ` ${this.sVarCalc}=${intSearch.toString()}`); let n = 0; let ok: boolean = false; let bAllowRestart = false; do { intSearch.setInterval(intSearch.intersect(intMax)); if (intSearch.length === 0) { this.debug(`searchTarget length= 0: ${this.sVarCalc}=${intSearch.toString()}`) break; } if (intSearch.hasTargetValue(rTarget)) { ok = true; break; } intSearch.growStep(2); intSearch.next(); if (bAllowRestart) { intSearch.checkDirection() } else if (n > this._startIntervalMaxSteps / 2) { bAllowRestart = true; intSearch.reInit();
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
} } while (n++ < this._startIntervalMaxSteps); this.debug(`searchTarget fin: Target ${this.analyticalSymbol}=${rTarget} dans` + ` ${this.sVarCalc}=${intSearch.toString()}`); return { ok, intSearch }; } /** * Détermine l'intervalle de recherche initial * @param rTarget valeur cible de la fonction * @param rInit valeur initiale de la variable */ private getStartInterval(rTarget: number, rInit: number): any { const prmDom: ParamDomain = this._paramX.domain; const min: number = prmDom.minValue; const max: number = prmDom.maxValue; if (prmDom.domain === ParamDomainValue.NOT_NULL) { if (rInit === 0) { rInit = 1e-8; } } const intMax: Interval = new Interval(min, max); try { intMax.checkValue(rInit); } catch (m) { return { ok: false, res: m }; } // sens de variation de la fonction const inc = this.isIncreasingFunction(rInit, intMax); let step: number = 0.001; if (BoolIdentity(this.CalculX(rInit) > rTarget, inc)) { step = -step; // par ex, la fonction est croissante et la valeur initiale // de la variable a une image par la fonction > valeur cible } // initialisation de l'intervalle de recherche const intSearch1: SearchInterval = new SearchInterval(this, rInit, step, this.DBG); // on cherche dans une première direction let a = this.searchTarget(rTarget, intSearch1, intMax); if (a.ok) { return a; } // il se peut que la fonction ne soit pas monotone et qu'il faille chercher dans l'autre direction const intSearch2 = new SearchInterval(this, rInit, -step, this.DBG); a = this.searchTarget(rTarget, intSearch2, intMax); if (a.ok) { return a; } // gestion de l'erreur // la valeur cible de la fonction est elle trouvable ? let res: Message; let errDomain = false; switch (prmDom.domain) { case ParamDomainValue.INTERVAL: const si: SearchInterval = new SearchInterval(this, intMax.min, intMax.max - intMax.min); errDomain = !si.hasTargetValue(rTarget); break; case ParamDomainValue.POS: case ParamDomainValue.POS_NULL: const y = this.CalculX(1e-8);
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
errDomain = inc && (rTarget < y); break; case ParamDomainValue.ANY: break; default: throw new Error("unsupported parameter domain value"); } if (errDomain) { res = new Message(MessageCode.ERROR_DICHO_INIT_DOMAIN); res.extraVar.targetSymbol = this.analyticalSymbol; // symbole de la variable calculée par la fonction res.extraVar.targetValue = rTarget; // valeur cible pour la fonction // intervalle de valeurs pour la variable d'entrée de la fonction res.extraVar.variableInterval = intMax.toString(); res.extraVar.variableSymbol = this._paramX.symbol; // symbole de la variable d'entrée de la fonction } else { // Fusion des infos des deux intervalles explorés const intFus: number[] = []; // Signification des indices : 0,1 : minmax target; 2, 3 : variables associées if (isNaN(intSearch1.targets.max) || isNaN(intSearch2.targets.max)) { // bug targets en NaN en prod mais pas en debug pas à pas en explorant les variables intSearch1.updateTargets(); intSearch2.updateTargets(); } if (intSearch1.targets.max > intSearch2.targets.max) { intFus[1] = intSearch1.targets.max; intFus[3] = intSearch1.getVal(intSearch1.targets.maxIndex); } else { intFus[1] = intSearch2.targets.max; intFus[3] = intSearch2.getVal(intSearch2.targets.maxIndex); } if (intSearch1.targets.min < intSearch2.targets.min) { intFus[0] = intSearch1.targets.min; intFus[2] = intSearch1.getVal(intSearch1.targets.minIndex); } else { intFus[0] = intSearch2.targets.min; intFus[2] = intSearch2.getVal(intSearch2.targets.minIndex); } if (intFus[1] < rTarget) { res = new Message(MessageCode.ERROR_DICHO_TARGET_TOO_HIGH); res.extraVar.extremeTarget = intFus[1]; res.extraVar.variableExtremeValue = intFus[3]; } else { res = new Message(MessageCode.ERROR_DICHO_TARGET_TOO_LOW); res.extraVar.extremeTarget = intFus[0]; res.extraVar.variableExtremeValue = intFus[2]; } res.extraVar.variableSymbol = this._paramX.symbol; // symbole de la variable de la fonction res.extraVar.targetSymbol = this.analyticalSymbol; // symbole de la variable calculée par la fonction res.extraVar.targetValue = rTarget; // valeur cible pour la fonction } return { ok: false, res }; } private CalcZero(x: number): number { this.vX = x; return this.Calcul() - this._target; } /** * Searches the interval from <tt>lowerLimit</tt> to <tt>upperLimit</tt> * for a root (i.e., zero) of the function <tt>func</tt> with respect to * its first argument using Brent's method root-finding algorithm. * * Translated from zeroin.c in http://www.netlib.org/c/brent.shar. * * Copyright (c) 2012 Borgar Thorsteinsson <borgar@borgar.net> * MIT License, http://www.opensource.org/licenses/mit-license.php
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
* * @param {function} func function for which the root is sought. * @param {number} lowerlimit the lower point of the interval to be searched. * @param {number} upperlimit the upper point of the interval to be searched. * @param {number} errorTol the desired accuracy (convergence tolerance). * @param {number} maxIter the maximum number of iterations. * @returns an estimate for the root within accuracy. * */ private uniroot<T>(func: (param: number) => number, thisArg: T, lowerLimit: number, upperLimit: number, errorTol: number = 0, maxIter?: number ): { found: boolean, value: number } { let a = lowerLimit; let b = upperLimit; let c = a; let fa = func.call(thisArg, a); let fb = func.call(thisArg, b); let fc = fa; let tolAct; // Actual tolerance let newStep; // Step at this iteration let prevStep; // Distance from the last but one to the last approximation let p; // Interpolation step is calculated in the form p/q; division is delayed until the last moment let q; if (maxIter === undefined) { maxIter = SessionSettings.maxIterations; } while (maxIter-- > 0) { prevStep = b - a; if (Math.abs(fc) < Math.abs(fb)) { // Swap data for b to be the best approximation a = b, b = c, c = a; fa = fb, fb = fc, fc = fa; } tolAct = 1e-15 * Math.abs(b) + errorTol / 2; newStep = (c - b) / 2; if ( Math.abs(newStep) <= tolAct && ( Math.abs(fb) < errorTol || Math.abs(fb - fc) < errorTol * 1E5 ) ) { // Acceptable approx. is found return { found: true, value: b }; } // Decide if the interpolation can be tried if (Math.abs(prevStep) >= tolAct && Math.abs(fa) > Math.abs(fb)) { // If prev_step was large enough and was in true direction, Interpolatiom may be tried let t1; let cb; let t2; cb = c - b; if (a === c) { // If we have only two distinct points linear interpolation can only be applied t1 = fb / fa; p = cb * t1; q = 1.0 - t1; } else { // Quadric inverse interpolation q = fa / fc, t1 = fb / fc, t2 = fb / fa; p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1)); q = (q - 1) * (t1 - 1) * (t2 - 1);
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458
} if (p > 0) { q = -q; // p was calculated with the opposite sign; make p positive } else { p = -p; // and assign possible minus to q } if (p < (0.75 * cb * q - Math.abs(tolAct * q) / 2) && p < Math.abs(prevStep * q / 2)) { // If (b + p / q) falls in [b,c] and isn't too large it is accepted newStep = p / q; } // If p/q is too large then the bissection procedure can reduce [b,c] range to more extent } if (Math.abs(newStep) < tolAct) { // Adjust the step to be not less than tolerance newStep = (newStep > 0) ? tolAct : -tolAct; } a = b, fa = fb; // Save the previous approx. b += newStep, fb = func.call(thisArg, b); // Do step to a new approxim. if ((fb > 0 && fc > 0) || (fb < 0 && fc < 0)) { c = a, fc = fa; // Adjust c for it to have a sign opposite to that of b } } // No acceptable approximation was found return { found: false, value: b }; } }