False-Position Method is one of Bisection Methods. If a real root is bounded by given two starting points, $$ xl $$ and $$ xu $$ of $$ f(x) = 0 $$, then we can approximate the solution by doing a linear interpolation between the points $$ [xl, f(xl)] $$ and $$ [xu, f(xu)] $$ to find the $$ xr $$ value.
Here is the source code of implementation of False-Position Method in C language.
1#define thres 0.001
2
3typedef struct {
4 int len;
5 float *coef;
6 int *exp;
7} poly;
Define a structure for polynomial.
len represents the number of non-zero terms.
coef represents an array of coefficients.
exp represents an array of exponents.
For example, if input is 4 1 4 -2 3 3 1 -5 0, it means $$ f(x) = x^4 - 2x^3 + 3x -5 $$.
And struct poly p has these data: len = 4, *coef = [1, -2, 3, -5], *exp = [4, 3, 1, 0]
1int main(void) {
2 ...
3 //read input polynomial from the file
4 fopen_s(&fi, "in.txt", "r");
5 fscanf_s(fi, "%f %f", &xl, &xu);
6 fscanf_s(fi, "%d", &n);
7
8 p.coef = malloc(n * sizeof(float));
9 p.exp = malloc(n * sizeof(int));
10
11 for (i = 0; i < n; i++) {
12 fscanf_s(fi, "%f %d", &p.coef[i], &p.exp[i]);
13 p.len += 1;
14 }
15
16 false_position(xl, xu, p);
17 free(p.coef);
18 free(p.exp);
19 fclose(fi);
20
21 return 0;
22}
The main function has a simple structure.
The procedure works in the function false_position().
1float calc_poly(poly p, float x) {
2 //return polynomial calculated value
3 float fx = 0.0;
4 int i;
5
6 for (i = 0; i < p.len; i++) {
7 fx += p.coef[i] * powf(x, (float)p.exp[i]);
8 }
9
10 return fx;
11}
This function computes the equation.
For example, when x = 1, the iteration processes $$ 1 _ (1)^4 + (-2)(1)^3 + 3 _ (1)^1 + (-5) $$ and the return value fx will be -3.
1int false_position(float xl, float xu, poly p) {
2 ...
3 do {
4 count++;
5
6 if (count > 1000) {
7 // prevent infinite loop
8 printf("False position: Failed\n");
9 return -1;
10 }
11
12 //procedure 1
13 fl = calc_poly(p, xl);
14 fu = calc_poly(p, xu);
15
16 if (fl * fu > 0) {
17 // when signs are same
18 printf("False position: Failed\n");
19 return -1;
20 }
21
22 //procedure2
23 xr = ((xl * fu) - (xu * fl)) / (fu - fl);
24 fr = calc_poly(p, xr);
25
26 //procedure3
27 if (fr < 0) {
28 xl = xr;
29 } else if (fr > 0) {
30 xu = xr;
31 } else if (fr == 0) {
32 //when find the root
33 break;
34 }
35
36 //procedure4
37 if (count == 1) {
38 xr_prev = xr;
39 } else if (count > 1) {
40 ea = (float)fabs((double)100 * ((xr - xr_prev) / xr));
41 xr_prev = xr;
42 }
43 } while (ea >= thres);
44
45 printf("False position: %f (%d iterations)\n", xr, count);
46 return 0;
47}
- Procedure 1 checks if there’s a root in the given range.
The root means a point where polynomial and x-axis meets, in other words, $$ f(x) = 0 $$.
If $$ f(x_l) \times f(x_u) > 0 $$, it means the signs are same(both positive or both negative), therefore no point of contact between
pand x-axis. - Procedure 2 calculates
xr, assumed root in this iteration, and the value $$ f(x_r) $$. - Procedure 3 determines the range for next iteration. This procedure reduces the range of possible root in the equation and increases the probability of finding the root.
- Procedure 4 calculates the error. When error is less then the threshold, the iteration terminates.
Comments