triangle_group_limit_set/exp_equation.c

#include "main.h"
#include "exp_equation.h"

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>

#define EPSILON 1e-9
#define LOOP(i) for(int i = 0; i < 3; i++)

struct solve_exp_plus_exp_params
{
	double alpha;
	double beta;
	double x;
};

static double solve_exp_plus_exp_f(double t, void *_params)
{
	struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;

//	return exp(params->alpha * t) + exp(params->beta * t) - params->x;
	return log(exp(params->alpha * t) + exp(params->beta * t)) - log(params->x);
}

static double solve_exp_plus_exp_df(double t, void *_params)
{
	struct solve_exp_plus_exp_params *params = (struct solve_exp_plus_exp_params*)_params;

//	return params->alpha * exp(params->alpha * t) + params->beta * exp(params->beta * t);
	return params->alpha + (params->beta - params->alpha) / (1 + exp((params->alpha-params->beta)*t));
}

static void solve_exp_plus_exp_fdf(double t, void *params, double *f, double *df)
{
	*f = solve_exp_plus_exp_f(t, params);
	*df = solve_exp_plus_exp_df(t, params);
}

// solve the equation exp(alpha t) + exp(beta t) = x for t
int solve_exp_plus_exp(double alpha, double beta, double x, double *t)
{
	if(alpha <= 0 && beta >= 0 || alpha >= 0 && beta <= 0) {
		double critical =
			pow(-beta/alpha, alpha/(alpha-beta)) +
			pow(-beta/alpha, beta/(alpha-beta));

		if(x < critical)
			return 0;
		else if (x < critical + EPSILON) {
			t[0] = log(-beta/alpha)/(alpha-beta);
			return 1;
		}

		// Newton this
		gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
		struct solve_exp_plus_exp_params params;
		params.alpha = alpha;
		params.beta = beta;
		params.x = x;
		gsl_function_fdf FDF;
		FDF.f = &solve_exp_plus_exp_f;
		FDF.df = &solve_exp_plus_exp_df;
		FDF.fdf = &solve_exp_plus_exp_fdf;
		FDF.params = (void *)&params;
		int status;
		double root, lastroot;

		for(int r = 0; r < 2; r++) {
			root = r == 0 ? log(x)/beta : log(x)/alpha;
			gsl_root_fdfsolver_set(solver, &FDF, root);
			for(int i = 0; i < 100; i++) {
				gsl_root_fdfsolver_iterate(solver);
				lastroot = root;
				root = gsl_root_fdfsolver_root(solver);
				status = gsl_root_test_delta(root, lastroot, 0, 1e-9);

				if(status == GSL_SUCCESS)
					break;

//				printf("iteration %d, root %f\n", i, root);
			}
			t[r] = root;
		}

		gsl_root_fdfsolver_free(solver);

		return 2;
	} else {
		// Newton with start value 0
		gsl_root_fdfsolver *solver = gsl_root_fdfsolver_alloc(gsl_root_fdfsolver_newton);
		struct solve_exp_plus_exp_params params;
		params.alpha = alpha;
		params.beta = beta;
		params.x = x;
		gsl_function_fdf FDF;
		FDF.f = &solve_exp_plus_exp_f;
		FDF.df = &solve_exp_plus_exp_df;
		FDF.fdf = &solve_exp_plus_exp_fdf;
		FDF.params = (void *)&params;
		int status;
		double root, lastroot;

		root = 0;
		gsl_root_fdfsolver_set(solver, &FDF, root);
		for(int i = 0; i < 100; i++) {
			gsl_root_fdfsolver_iterate(solver);
			lastroot = root;
			root = gsl_root_fdfsolver_root(solver);
			status = gsl_root_test_delta(root, lastroot, 0, 1e-9);

//			printf("iteration %d, root %f\n", i, root);

			if(status == GSL_SUCCESS)
				break;
		}
		t[0] = root;

		gsl_root_fdfsolver_free(solver);

		return 1;
	}
}

// solve the equation x1 exp(a1 t) + x2 exp(a2 t) + x3 exp(a3 t) = 0
int solve_linear_exp(vector_t a, vector_t x, double *t)
{
	if(x.x[0] > 0 && x.x[1] > 0 && x.x[2] > 0 ||
	   x.x[0] < 0 && x.x[1] < 0 && x.x[2] < 0)
		return 0;

	// ensure that y[0] < 0 and y[1], y[2] > 0
	int j;
	vector_t y, b;
	for(j = 0; j < 3; j++)
		if(x.x[(j+1)%3] * x.x[(j+2)%3] > 0)
			break;
	LOOP(i) y.x[i] = x.x[(i+j)%3];
	if(y.x[0] > 0)
		LOOP(i) y.x[i] *= -1;
	LOOP(i) b.x[i] = a.x[(i+j)%3];

	double T = (log(y.x[1]) - log(y.x[2])) / (b.x[2] - b.x[1]);
	double rhs = - y.x[0] *
		pow(y.x[1], (b.x[0]-b.x[2])/(b.x[2]-b.x[1])) *
		pow(y.x[2], (b.x[0]-b.x[1])/(b.x[1]-b.x[2]));

	int n = solve_exp_plus_exp(b.x[1] - b.x[0], b.x[2] - b.x[0], rhs, t);

	for(int i = 0; i < n; i++)
		t[i] += T;

	return n;
}

/*
int main(int argc, char *argv[])
{

//	int n = solve_exp_plus_exp(atof(argv[1]), atof(argv[2]), atof(argv[3]), result);

	double result[2];
	vector_t a, x;
	a.x[0] = atof(argv[1]);
	a.x[1] = atof(argv[2]);
	a.x[2] = atof(argv[3]);
	x.x[0] = atof(argv[4]);
	x.x[1] = atof(argv[5]);
	x.x[2] = atof(argv[6]);
	int n = solve_linear_exp(a, x, result);

	if(n == 0)
		printf("0 results found\n");
	else if(n == 1)
		printf("1 result found: %.9f\n", result[0]);
	else if(n == 2)
		printf("2 results found: %.9f and %.9f\n", result[0], result[1]);
	else
		printf("%d results found\n", n);
	return 0;
}
*/