[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] One dimensional root of an equation
From: |
Rajil Saraswat |
Subject: |
[Help-gsl] One dimensional root of an equation |
Date: |
Sat, 30 Oct 2004 18:53:32 +0100 |
User-agent: |
KMail/1.7.1 |
Hi,
I am trying to solve a quadratic equation given as follows:
y=441682.9251*(-1.842052048*x+0.3941686334e-2)/((-x+0.2e-3)*(-525.00*x+2.716111111))-242538.2570*(1.953506658*x+0.7436592362e-3)/(x*(525.00*x+1.361111111));
I tried maple which gave three roots to me {x = -0.001365774428}, {x =
0.00003685022686}, {x = 0.002880307005} out of which the root i am looking
for is x = 0.00003685022686.
So far so good, next i tried gsl with the brent algorithm as i wanted to
write a C program for interfacing purposes. However the program failed with
the error:
gsl: brent.c:57: ERROR: function not continuous
Default GSL error handler invoked.
Aborted
Is there any way (gsl or otherwise) i can solve this equation? Here is the
gsl code for reference.
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
struct quadratic_params
{
double a, b, c;
};
double
quadratic (double x, void *params)
{
struct quadratic_params *p
= (struct quadratic_params *) params;
const double y0=
441682.9251*(-1.842052048*x+0.3941686334e-2)/((-x+0.2e-3)*(-525.00*x+2.716111111))-242538.2570*(1.953506658*x+0.7436592362e-3)/(x*(525.00*x+1.361111111));
return y0;
}
int
main (void)
{
int status;
int iter = 0, max_iter = 100;
const gsl_root_fsolver_type *T;
gsl_root_fsolver *s;
double r = 0, r_expected = 0.3685022686e-4;
double x_lo = 0.0, x_hi =0.0001;
gsl_function F;
struct quadratic_params params = {1.0, 0.0, -5.0};
F.function = &quadratic;
F.params = ¶ms;
T = gsl_root_fsolver_brent;
s = gsl_root_fsolver_alloc (T);
gsl_root_fsolver_set (s, &F, x_lo, x_hi);
printf ("using %s method\n",
gsl_root_fsolver_name (s));
printf ("%5s [%9s, %9s] %9s %10s %9s\n",
"iter", "lower", "upper", "root",
"err", "err(est)");
do
{
iter++;
status = gsl_root_fsolver_iterate (s);
r = gsl_root_fsolver_root (s);
x_lo = gsl_root_fsolver_x_lower (s);
x_hi = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval (x_lo, x_hi,
0, 0.001);
if (status == GSL_SUCCESS)
printf ("Converged:\n");
printf ("%5d [%.7f, %.7f] %.7f %+.7f %.7f\n",
iter, x_lo, x_hi,
r, r - r_expected,
x_hi - x_lo);
}
while (status == GSL_CONTINUE && iter < max_iter);
return status;
}
- [Help-gsl] One dimensional root of an equation,
Rajil Saraswat <=