I have been trying to code the ujevic method for solving non linear equations. but upon compiling the code with given data, it return the root as - nan. The actual root should be 1 and no. of iterations 75.
#include <stdio.h>
#include <math.h>
#define f(x) ((x*x*x)-(3*x) 2)
#define g(x) ((3*x*x)-3)
void main() {
float a, b, x0, xk, z0, xt, y;
int i = 0;
printf("\nEnter The Value Of a: ");
scanf("%f", &a);
printf("\nEnter The Value Of b: ");
scanf("%f", &b);
printf("\nEnter The Value Of x0: ");
scanf("%f", &x0);
printf("\nEnter The Value Of Alpha: ");
scanf("%f", &y);
do {
xt = x0;
z0 = x0 - y * (f(x0) / g(x0));
xk = x0 4 * (z0 - x0) * (f(x0) / ((3 * f(x0)) - (2 * f(z0))));
x0 = xk;
i ;
} while (fabs(xk - xt >= 0.000000000000000001));
printf("The value of root is: %f", xk);
printf("\n the no of iterations is %d", i);
}
Output:
Enter The Value Of a: 0.5
Enter The Value Of b: 1.2
Enter The Value Of x0: 0.5
Enter The Value Of Alpha: 0.1
The value of root is: -nan
the no of iterations is 43
CodePudding user response：
When x0 reaches the value 0.999802410602569580078125, f(x0) is calculated as zero. In f(x0), x0*x0*x0 is calculated as 0.999407351016998291015625, and 3*x0 is 2.999407291412353515625. Subtracting these would give −1.999999940395355224609375 mathematically, but the float format in your C implementation (likely IEEE-754 binary32) cannot represent that number. The closest numbers it can represent are −1.99999988079071044921875 and −2. Of those, −2 is closer to −1.999999940395355224609375 than −1.99999988079071044921875 is, so, when x0*x0*x0 - 3*x0 is calculated using float arithmetic, the result is −2. Then adding 2, the final step in the f macro, yields zero.
The calculation of f(z0) proceeds similarly, and then f(x0)/((3*f(x0))-(2*f(z0))) evaluates as 0 / (0-0), which yields a NaN.
You have two mathematical problems in the code, aside from various other issues. One is that your loop approaches a point where f(x0)/((3*f(x0))-(2*f(z0))) will be mathematically zero. At some point, the calculations in that loop will fail. I have not looked at the Ujevic method and do not know what the intent is here. At the very least, once ((3*f(x0))-(2*f(z0))) evaluates as zero, you cannot continue the loop anymore.
Second, fabs(xk-xt) >=0.000000000000000001), which you apparently intended to write instead of fabs(xk-xt>=0.000000000000000001) is futile with float variables using the common format for float whenb xk is near 1. In this format, the smallest different between two numbers just under 1 is 0.000000059604644775390625, which is much larger than 0.000000000000000001. It is therefore impossible for fabs(xk-xt) to be smaller than 0.000000000000000001 unless it is zero, meaning xk and xt are equal. Even if the variables used the format commonly used for double (IEEE-754 binary64), 0.000000000000000001 is smaller than the smallest difference representable in numbers just under 1. You need a larger bound in this test, or you need to ensure the loop will always converge to identical xk and xt.
CodePudding user response：
There are multiple problems in your code:
the return type of
mainisint, notvoid.you should use type
doubleinstead offloat: the limited precision of typefloatcauses errors with your computation.the test
while(fabs(xk - xt >= 0.000000000000000001))is incorrect. You probably mean:while (fabs(xk - xt) >= 0.000000000000000001)the argument
xin the expansion of macrosf(x)andg(x)should be parenthesized to avoid operator precedence issues when passing non trivial expressions as macro arguments.
Here is a modified version:
#include <stdio.h>
#include <math.h>
#define f(x) ((x)*(x)*(x)-3*(x) 2)
#define g(x) (3*(x)*(x)-3)
int main() {
double a, b, x0, xk, z0, xt, y;
int i = 0;
printf("\nEnter The Value Of a: ");
if (scanf("%lf", &a) != 1)
return 1;
printf("\nEnter The Value Of b: ");
if (scanf("%lf", &b) != 1)
return 1;
printf("\nEnter The Value Of x0: ");
if (scanf("%lf", &x0) != 1)
return 1;
printf("\nEnter The Value Of Alpha: ");
if (scanf("%lf", &y) != 1)
return 1;
do {
xt = x0;
z0 = x0 - y * (f(x0) / g(x0));
xk = x0 4 * (z0 - x0) * (f(x0) / ((3 * f(x0)) - (2 * f(z0))));
x0 = xk;
i ;
} while (fabs(xk - xt) >= 1.0e-18);
printf("The value of root is: %f\n", xk);
printf("the number of iterations is %d\n", i);
return 0;
}
The output is still nan:
Enter The Value Of a: 0.5
Enter The Value Of b: 1.2
Enter The Value Of x0: 0.5
Enter The Value Of Alpha: 0.1
The value of root is: nan
the number of iterations is 101
Using the more precise type long double, there are more iterations, but still the same nan output:
Enter The Value Of a: 0.5
Enter The Value Of b: 1.2
Enter The Value Of x0: 0.5
Enter The Value Of Alpha: 0.1
The value of root is: nan
the number of iterations is 119
Your algorithm fails because (3 * f(x0) - 2 * f(z0)) becomes zero before the xt and xk get closer than 1.0e-18.
Adding a sanity test if (3 * f(x0) - 2 * f(z0) == 0) break; and initializing xk = x0; before the loop gives this output:
Enter The Value Of a: 0.5
Enter The Value Of b: 1.2
Enter The Value Of x0: 0.5
Enter The Value Of Alpha: 0.1
The value of root is: 1.000000
the number of iterations is 100
Note that the convergence test is incorrect: fabs(xk - xt) >= 1.0e-18 tests the absolute difference, not the relative difference. You should instead use something like:
while (fabs(xk - xt) / fmax(fabs(xk), fabs(xt)) >= 1.0e-8);
Here is an updated version that does converge in fewer than 100 iterations, ie: before the sanity test kicks:
#include <stdio.h>
#include <math.h>
#define f(x) ((x)*(x)*(x)-3*(x) 2)
#define g(x) (3*(x)*(x)-3)
int main() {
double a, b, x0, xk, z0, xt, y;
int i = 0;
printf("\nEnter The Value Of a: ");
if (scanf("%lf", &a) != 1)
return 1;
printf("\nEnter The Value Of b: ");
if (scanf("%lf", &b) != 1)
return 1;
printf("\nEnter The Value Of x0: ");
if (scanf("%lf", &x0) != 1)
return 1;
printf("\nEnter The Value Of Alpha: ");
if (scanf("%lf", &y) != 1)
return 1;
xk = x0;
do {
xt = x0;
z0 = x0 - y * (f(x0) / g(x0));
if (3 * f(x0) - 2 * f(z0) == 0)
break;
xk = x0 4 * (z0 - x0) * (f(x0) / ((3 * f(x0)) - (2 * f(z0))));
x0 = xk;
i ;
} while (fabs(xk - xt) / fmax(fabs(xk), fabs(xt)) >= 1.0e-8);
printf("The value of root is: %f\n", xk);
printf("the number of iterations is %d\n", i);
return 0;
}
Output:
Enter The Value Of a: 0.5
Enter The Value Of b: 1.2
Enter The Value Of x0: 0.5
Enter The Value Of Alpha: 0.1
The value of root is: 1.000000
the number of iterations is 88