I'm trying to built a function which i can give as argument a mathematical expression such as "x**3 bx-c" and trough the newton method return me the roots of the given function.
I was thinking in something like this following example, where 'r0' is an initial random number non-zero and dfunc is the derivative of the function. I don't know how i'll evaluate the derivative, but i want to know first if is possible what i want to do with the expression as an argument.
double * roots (double r0, <type> func) {
double *r = alloc(100, sizeof(double));
r[0] = r0;
for (int i = 1; i <= 100; i ) {
r[i]=r[i-1]-func(r[i-1])/dfunc(r[i-1]);
};
return r;
};
CodePudding user response：
Given OP's function is like double function(double x), pass the address of that function to roots().
Create a derivative routine that calls func(x) from 2 nearby values to determine the slope.
OP's for loop only iterates to find 1 root. Additional code needed to find more.
Also no need for to allocate an array of double when looking for one root. Enough to keep track of current x and previous x.
extern func(double func(double), x);
double root(double r0, double func(double x)) {
double x =r0;
for (int i = 1; i < 100; i ) {
double better_x = x - func(x)/dfunc(func, x);
x = better_x;
};
return x;
}
// Usage
double y = root(0.0, function);
OP's application of Netwon's method deserves detection of a zero slope, excessive iterations, better stopping criteria than a simple iteration count, etc.
Further, for a general double function(double x) (and not just a polynomial one) to find the N 1 roots, additional code needed.
IMO, OP show instead focus on polynomial roots only and create a function of the form:
complex double *polynomial_roots(size_t n, const double coefficients[n 1]);
See Bairstow's method.
IAC, OP's goal multi-root goal is non-trivial.
CodePudding user response：
Others have suggested using function pointes, and I agree that this is probably the correct solution, but you did ask for a way to provide expressions to a function. You can do this, of course, but you have to implement expressions yourself.
You can make a structure for the various kinds of sub-expressions you want, and it could look like this:
struct expr
{
int refcount; // GC; exprs share sub-expressions
enum
{
CONS, VAR, MINUS, /* unary minus */
ADD, SUB, /* binary subtraction */
MUL, /* more... */
} type;
union {
struct { double x; } cons;
struct { int idx; /* index in var vector */ } var;
struct { struct expr *x; } unary;
struct { struct expr *x, *y; } binary;
};
};
I've just stuck to some easy ones here.
In the derivative code below, I need to keep expressions in the multiplication rule (and you would need to do the same with the chain rule if you implement that), so I handle memory management with reference counters.
struct expr *incref(struct expr *expr)
{
expr->refcount ;
return expr;
}
void free_expr(struct expr *);
struct expr *decref(struct expr *expr)
{
expr->refcount--;
if (expr->refcount > 1)
return expr;
free_expr(expr);
return 0;
}
I return a pointer to the object I modify because I will need that for incref. I don't use it for decref, but I preferred to keep the interface for the two functions consistent.
There is nothing much to creating instances of expressions; you just allocate and set the values in them.
struct expr *new_const(double x)
{
struct expr *expr = malloc(sizeof *expr);
assert(expr); // handle alloc errors
*expr = (struct expr){.refcount = 1, .type = CONS, .cons = {x}};
return expr;
}
struct expr *new_var(int idx)
{
struct expr *expr = malloc(sizeof *expr);
assert(expr); // handle alloc errors
*expr = (struct expr){.refcount = 1, .type = VAR, .var = {idx}};
return expr;
}
struct expr *new_minus(struct expr *x)
{
struct expr *expr = malloc(sizeof *expr);
assert(expr); // handle alloc errors
*expr = (struct expr){.refcount = 1, .type = MINUS, .unary = {x}};
return expr;
}
struct expr *new_add(struct expr *x, struct expr *y)
{
struct expr *expr = malloc(sizeof *expr);
assert(expr); // handle alloc errors
*expr = (struct expr){.refcount = 1, .type = ADD, .binary = {x, y}};
return expr;
}
struct expr *new_sub(struct expr *x, struct expr *y)
{
struct expr *expr = malloc(sizeof *expr);
assert(expr); // handle alloc errors
*expr = (struct expr){.refcount = 1, .type = SUB, .binary = {x, y}};
return expr;
}
struct expr *new_mul(struct expr *x, struct expr *y)
{
struct expr *expr = malloc(sizeof *expr);
assert(expr); // handle alloc errors
*expr = (struct expr){.refcount = 1, .type = MUL, .binary = {x, y}};
return expr;
}
For functions that work on them, you can use a switch:
void print_expr(struct expr *expr)
{
switch (expr->type)
{
case CONS:
printf("%f", expr->cons.x);
break;
case VAR:
printf("x%d", expr->var.idx);
break;
case MINUS:
printf("-(");
print_expr(expr->unary.x);
printf(")");
break;
case ADD:
printf("(");
print_expr(expr->binary.x);
printf(" ");
print_expr(expr->binary.y);
printf(")");
break;
case SUB:
printf("(");
print_expr(expr->binary.x);
printf(" - ");
print_expr(expr->binary.y);
printf(")");
break;
case MUL:
printf("(");
print_expr(expr->binary.x);
printf(" * ");
print_expr(expr->binary.y);
printf(")");
break;
/* more here ... */
}
}
With function pointers you can implement polymorphism to avoid the switching, but that is probably overkill here.
If you want to evaluate an expression, you need to provide values for the variables, and an easy solution is to just use an array. That is why the variables above contain an index rather than a string or similar. The index is interpreted as the offset into the array of variable values.
double eval_expr(struct expr *expr, double *x)
{
switch (expr->type)
{
case CONS:
return expr->cons.x;
case VAR:
return x[expr->var.idx];
case MINUS:
return -eval_expr(expr->unary.x, x);
case ADD:
return eval_expr(expr->binary.x, x) eval_expr(expr->binary.y, x);
case SUB:
return eval_expr(expr->binary.x, x) - eval_expr(expr->binary.y, x);
case MUL:
return eval_expr(expr->binary.x, x) * eval_expr(expr->binary.y, x);
/* more here ... */
}
}
To calculate the derivative, you implement transformation rules according to the derivation rules (which is why I kept to simple expressions here; I didn't want to write code for the chain rule or division and such; I am a very lazy man).
struct expr *derivative(struct expr *expr, int var)
{
switch (expr->type)
{
case CONS:
return new_const(0);
case VAR:
return (expr->var.idx == var) ? new_const(1) : new_const(0);
case MINUS:
return new_minus(derivative(expr->unary.x, var));
case ADD:
return new_add(derivative(expr->binary.x, var), derivative(expr->binary.y, var));
case SUB:
return new_sub(derivative(expr->binary.x, var), derivative(expr->binary.y, var));
case MUL:
return new_add(new_mul(incref(expr->binary.x), derivative(expr->binary.y, var)),
new_mul(derivative(expr->binary.x, var), incref(expr->binary.y)));
/* more here ... */
}
}
I use incref in the multiplication rule, so I can reuse the incoming expressions and still free memory safely. Freeing an expression is quite simple then:
void free_expr(struct expr *expr)
{
switch (expr->type)
{
case CONS:
case VAR:
break;
case MINUS:
decref(expr->unary.x);
break;
case ADD:
case SUB:
case MUL:
decref(expr->binary.x);
decref(expr->binary.y);
break;
/* more here ... */
}
free(expr);
}
You can use these expressions like this:
int main(void)
{
struct expr *expr =
new_mul(
new_add(
new_const(3.5),
new_minus(new_var(0))),
new_sub(new_const(1.4),
new_var(1)));
print_expr(expr);
putchar('\n');
double x[] = {2.0, 42.0}; // variables
printf("value: %f\n", eval_expr(expr, x));
struct expr *deriv = derivative(expr, 0);
printf("derivative with respect to x0:\n");
print_expr(deriv);
putchar('\n');
printf("derivative value: %f\n", eval_expr(deriv, x));
//free_expr(deriv);
deriv = derivative(expr, 1);
printf("derivative with respect to x1:\n");
print_expr(deriv);
putchar('\n');
printf("derivative value: %f\n", eval_expr(deriv, x));
free_expr(deriv);
free_expr(expr);
return 0;
}
There are most likely errors in this code, I just whipped it up quickly, but it should give an idea about how you could work with expressions.