commercial/solvers.cpp

// $Id: solvers.cpp,v 1.3 2008/02/01 06:24:29 d3g637 Exp $
// solver.cpp: general equation solver for compound exponential decays
//

#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <time.h>

#include "solvers.h"

#ifndef WIN32
#define _finite isfinite
#endif

#define SGNP(X,P) ((X)<-(P)?-1:((X)>(P)?+1:0))

int dual_decay_solve(double *ans, double prec, double start, double end, int f, double a, double n, double b, double m, double c)
{
#define EVAL(f,t,a,n,b,m,c) (a*pow(n,f)*exp(n*t) + b*pow(m,f)*exp(m*t) + c)

    // check differential for sign change over range
    double df_start = EVAL(f+1,start,a,n,b,m,c);
    double df_end = EVAL(f+1,end,a,n,b,m,c);

    if(!_finite(df_start) || !_finite(df_end))
        return -1;

    int sa = SGNP(df_start,prec);
    int sb = SGNP(df_end,prec);

    // if sign of first differential changed over range
    if (sa*sb==-1)
    {
        // min/max in start-end range implies multiple solutions
        double t;
        if (!dual_decay_solve(&t, prec, start, end, a*pow(n,f+1),n,b*pow(m,f+1),m,0))

            // error: should have a solution for the inflexion point
            return -1;
        
        // look for solutions right and left of first min/max, keep left result
        int n_right = dual_decay_solve(ans, prec, t, end, f, a,n,b,m,c);
        if (n_right<0)
            return -1; // error; right solution failed
        
        int n_left = dual_decay_solve(ans, prec, start, t, f, a,n,b,m,c);
        if (n_left<0)
            return -1; // error: left solution failed

        return n_left+n_right;
    }
    return dual_decay_solve(ans, prec, start, end, a,n,b,m,c);
#undef EVAL
}

int dual_decay_solve(double *ans, double prec, double start, double end, double a, double n, double b, double m, double c)
{
#define EVAL(t,a,n,b,m,c) (a*exp(n*t) + b*exp(m*t) + c)
    // check for sign change over range
    double f_start = EVAL(start,a,n,b,m,c);
    double f_end = EVAL(end,a,n,b,m,c);

    if (!_finite(f_start) ||!_finite(f_end))
        return -1;

    int sa = SGNP(f_start,prec);
    int sb = SGNP(f_end,prec);

    // if sign hasn't change in range
    if (sa*sb==1)

        // no solution possible
        return 0;

    // check for boundary solutions
    else if (sa==0)
        *ans = start;
    else if (sb==0)
        *ans = end;
    else // opposite signs
    {
        // only one solution exists -- use Newton's method to find it
        double t=(start+end)/2;
        int result = dual_decay_solve(ans,prec,start,t,a,n,b,m,c);
        if (result==0)
            result = dual_decay_solve(ans,prec,t,end,a,n,b,m,c);
        return result;
    }   

    return 1;
}

#undef SGNP

---

GridLAB-D^TM Version 1.0

An open-source project initiated by the US Department of Energy