Public feed
#722492
#include <stdio.h>
#include <math.h>
#include "Problem6.c" // Includes the helper functions for the actual Problem
#include "TestCase1.c" // Includes the helper functions for square root of 2
// A macro to ease printing the roots for the results
#define PRINT_ROOT(a) printf("root: %.18Le\n", (long double)a);
// Epsilon as 10^-18
#define eps (long double)0.000000000000000001
// Function Prototypes, so that main could be moved up, works better for
// Commenting all the code.
long double Bisection(long double (*f)(long double), long double a, long double b);
short sign(long double x);
long double Bisection(long double (*f)(long double), long double a, long double b);
long double Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0);
long double Secant(long double (*f)(long double), long double a, long double b);
long double Multi_Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0, int p);
// Main function, calls the necessary functions for the project
// Liberal use of function pointers follows. Quick reference:
// &funtion_name returns a pointer to a functions
// (returntype)(*var)(paramters) is the type for a function pointer saved in var
// (*var)(arguments) calls a function pointer stored in var.
// -----------------------------------------------------------------------------
// Implementations of the three root finding methods follow this function
int main() {
long double root;
printf("Numerical Analysis Project 2, Problem 6 Approximator\n");
printf("Using epsilon = %.11Le\n\n", eps);
// Calls the test case, from TestCase1.c
printf("Test Case, 0=x^2-2, guess between 0 and 3, or 1.5 for Newton's method\n");
root = Bisection(&x_squared_minus_2, 0.0, 3.0);
PRINT_ROOT(root);
root = Newton(&x_squared_minus_2, &x_squared_minus_2_prime, 1.5);
PRINT_ROOT(root);
root = Secant(&x_squared_minus_2, 0.0, 3.0);
PRINT_ROOT(root);
printf("\n");
// Solves the problem, Problem 6 from the sheet, Electrical Engineering
printf("Problem 6, unintentional guesses of 0 and 3, or 1.5 for Newton's method\n");
root = Bisection(&Problem6, 0.0, 3.0);
PRINT_ROOT(root);
root = Newton(&Problem6, &Problem6_prime, 1.5);
PRINT_ROOT(root);
root = Multi_Newton(&Problem6, &Problem6_prime, 1.5, 2);
PRINT_ROOT(root);
root = Secant(&Problem6, 0.0, 3.0);
PRINT_ROOT(root);
printf("\n");
printf("Problem 6, actual guesses of 1*10^10 and 2*10^10, or 1.5*10^10 for Newton\n");
root = Bisection(&Problem6, 1.0*powl(10.0, 10.0), 2.0*powl(10.0, 10.0));
PRINT_ROOT(root);
root = Newton(&Problem6, &Problem6_prime, 1.5*powl(10.0, 10.0));
PRINT_ROOT(root);
root = Multi_Newton(&Problem6, &Problem6_prime, 1.5*powl(10.0, 10.0), 2);
PRINT_ROOT(root);
root = Secant(&Problem6, 1.0*powl(10.0, 10.0), 2.0*powl(10.0, 10.0));
PRINT_ROOT(root);
return 0;
}
// Utility function to determine the sign of a long double
// returns 1 for positive or 0 and -1 for < 0
short sign(long double x) {
if(x < 0.0)
return -1;
return 1;
}
// Implementation of the Bisection method
long double Bisection(long double (*f)(long double), long double a, long double b) {
long double c; // Center value
int i = 1; // Counter value to measure the speed of convergence
// Loop that runs forever, incrementing i on each iteration
for(;;i++) {
// Calculate the new center value for this iteration
c = (a + b) / 2.0;
// If the difference between b and c is < epsilon we found the root
if((b-c) < eps) {
printf("Bisection method in %d iterations\n", i);
return c;
}
// if the root is between b and c
if( sign(((*f)(c)) * sign((*f)(b))) < 0)
// set a to be c
a = c;
else
// set b to be c
b = c;
// Next iteration
}
}
// Implementation of Newton's Method
long double Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0) {
int i = 1; // Counter value to prevent Newton from running
// too many iterations
long double x_1 = 0.0; // next x value
long double fprime_of_x_0; // Derivative of f at the current x value
// Loop that runs forever, incrementing i each time
for(;;i++) {
// Calculate the derivative of the current x value
fprime_of_x_0 = (*fprime)(x_0);
// If the derivative of f is 0 at the current x value, Newton fails
if(fabsl(fprime_of_x_0) < eps) {
printf("WARN: Newton fails, f'(x_0) = 0\n");
return 0.0;
}
// Calculate the next x value
x_1 = x_0 - ((*f)(x_0) / fprime_of_x_0);
// If the difference between the current and next x values is very small
// We have found the root.
if(fabsl(x_1 - x_0) < eps) {
printf("Newton's method in %d iterations\n", i);
return x_1;
}
// Counter case to ensure Newton's Method doesn't run forever
if(i > 150) { // Default 50 iterations before failure, bumped up
printf("WARN: Newton took too long to converge, quitting\n");
return x_1;
}
//Set the current x value equal to the initial for the next iteration.
x_0 = x_1;
}
}
// Implementation of Newton's Method for Multiple Roots
long double Multi_Newton(long double (*f)(long double), long double (*fprime)(long double), long double x_0, int p) {
int i = 1; // Counter value to prevent Newton from running
// too many iterations
long double x_1 = 0.0; // next x value
long double fprime_of_x_0; // Derivative of f at the current x value
// Loop that runs forever, incrementing i each time
for(;;i++) {
// Calculate the derivative of the current x value
fprime_of_x_0 = (*fprime)(x_0);
// If the derivative of f is 0 at the current x value, Newton fails
if(fabsl(fprime_of_x_0) < eps) {
printf("WARN: Newton fails, f'(x_0) = 0\n");
return 0.0;
}
// Calculate the next x value
x_1 = x_0 - ((long double)p*(*f)(x_0) / fprime_of_x_0);
// If the difference between the current and next x values is very small
// We have found the root.
if(fabsl(x_1 - x_0) < eps) {
printf("Newton's method for p=%d in %d iterations\n", p, i);
return x_1;
}
// Counter case to ensure Newton's Method doesn't run forever
if(i > 150) { // Default 50 iterations before failure, bumped up
printf("WARN: Newton took too long to converge, quitting\n");
return x_1;
}
//Set the current x value equal to the initial for the next iteration.
x_0 = x_1;
}
}
// Implementation of the Secant Method
long double Secant(long double (*f)(long double), long double a, long double b) {
long double c; // Center value
int i = 1; // Counter value to measure the speed of convergence and
// prevent secant method from running forever, no guarantee of
// convergence.
// Loop that runs forever, incrementing i each time
for(;;i++) {
// Calculate the center value by the secant method
c = b - ((long double)(b - a)*(*f)(b)) / ((*f)(b) - (*f)(a));
// If the difference between the center and the upper bound is small
// We have found the root
if(fabsl(c - b) < eps) {
printf("Secant method in %d iterations\n", i);
return c;
}
// Counter case to ensure the Secant Method doesn't run forever
if(i > 50) { // Default 50 iterations before failure, bumped up
printf("WARN: Secant took too long to converge, quitting\n");
return c;
}
// Set the new bounds for the next iteration
a = b;
b = c;
}
}