Write a program that implements the secant method for root finding. Your program

Question

Write a program that implements the secant method for root finding. Your program is a variant of the Newton-Raphson method, so you should be able to modify the Chapter 6_11 code from Etter & Ingber.

Here is the code to modify:

/*----------------------------------------------------------*/
/* Program chapter6_11 */
/* */
/* This program finds the real roots of a cubic polynomial */
/* using the Newton-Raphson method. */
#include<iostream> //Required for cin, cout
#include<cmath> //Required for pow()
using namespace std;
int main()
{
// Declare objects.
int iterations(0);
double a0, a1, a2, a3, x, p, dp, tol;
// Get user input.
cout << "Enter coefficients a0, a1, a2, a3 ";
cin >> a0 >> a1 >> a2 >> a3;
cout << "Enter initial guess for root ";
cin >> x;
// Evaluate p at initial guess.
p = a0*pow(x,3) + a1*x*x + a2*x + a3;
// Determine tolerance.
tol = fabs(p);
while(tol > 0.001 && iterations < 100)
{
// Calculate the derivative.
dp = 3*a0*x*x + 2*a1*x + a2;
// Calculate next estimated root.
x = x - p/dp;
// Evaluate p at estimated root.
p = a0*x*x*x + a1*x*x + a2*x + a3;
tol = fabs(p);
iterations++;
}
if(tol < 0.001)
{
cout << "Root is " << x << endl
<< iterations << " iterations ";
}
else
cout << "Did not converge after 100 iterations ";
return 0;
}

Use your programs for the following physics problem:

Modern Physics: The black body radiation curve has a peak that can be determined by maximizing the Planck radiation formula in frequency form with respect to frequency. The Planck radiation formula is

If we take the derivative of the equation above and set it equal to zero and then substitute x = (hn/kT) we get an expression that reduces to

Solve for the roots of this equation to determine the constant for Wien’s law. Use the bisection method and the Newton Raphson method. You will need to replace the function in each program with equation 1 above.

From your value of x determine the constant in Wien’s law:

Add a few lines of code to each program to (a) calculate the value of the constant and (b) compare your results to the standard constant (i.e. take a % error) from a your textbook.

Explanation / Answer

Secant method is defined as:

xn = xn-1 - f(xn-1)*((xn-1-xn-2)/(f(xn-1) - f(xn-2)))

Implementation of secant method for root finding of cubic polynomial is below. Here x2, x1, p2, p1 are for xn-1, xn-2, f(xn-1), f(xn-2), respectively.

d is f(xn-1)*((xn-1-xn-2)/(f(xn-1) - f(xn-2))).

//////// Code starts here ///////////

#include <iostream>
#include<cmath>

using namespace std;
int main() {
// Declare objects.
int iterations(0);
double a0, a1, a2, a3, x, x1, x2, p1, p2, d, tol;
// Get user input.
cout << "Enter coefficients a0, a1, a2, a3 ";
cin >> a0 >> a1 >> a2 >> a3;
cout << "Enter initial two guesses for root ";
cin >> x1 >> x2;
// Evaluate p at initial guesses.
p1 = a0*pow(x1,3) + a1*x1*x1 + a2*x1 + a3;
p2 = a0*pow(x2,3) + a1*x2*x2 + a2*x2 + a3;
// Determine tolerance.
tol = fabs(p2);
while(tol > 0.001 && iterations < 100)
{
d = p2*((x2-x1)/(p2-p1));
// Calculate next estimated root.
x = x2 - d;
x1 = x2;
x2 = x;
// Evaluate p at estimated root.
p1 = a0*pow(x1,3) + a1*x1*x1 + a2*x1 + a3;
p2 = a0*pow(x2,3) + a1*x2*x2 + a2*x2 + a3;
tol = fabs(p2);
iterations++;
}
if(tol < 0.001)
{
cout << "Root is " << x << endl
<< iterations << " iterations ";
}
else
cout << "Did not converge after 100 iterations ";
return 0;
}

//////// Code ends here ///////////

Sample Output:

Enter coefficients a0, a1, a2, a3
4 5 3 6

Enter initial two guesses for root
3 4
Root is -1.44774

13 iterations

Physics part:

Code to find x using Newton-Raphson method:

Note: This program is not taking any input at present. All the constants values are specified in.

//////// Code starts here ///////////

#include<iostream>
#include<cmath>
using namespace std;
int main()
{
// Declare objects.
int iterations(0);
double h, c, k, x, p, dp, tol;
// Get user input.
//cout << "Enter the values of constants h, c, k ";
//cin >> h >> c >> k >> a3;
//cout << "Enter initial guess for root ";
//cin >> x;
// Evaluate p at initial guess.
h = 6.626176*pow(10, -34);
c = 299792458;
k = 1.38064852*pow(10, -23);
x = 20;
  
//x+3e^{-x} -3=0
p = 3*exp(-x) + x - 3;
// Determine tolerance.
tol = fabs(p);
while(tol > 0.00001 && iterations < 100)
{
// Calculate the derivative.
dp = (-3)*exp(-x) + 1;
// Calculate next estimated root.
x = x - p/dp;
// Evaluate p at estimated root.
p = 3*exp(-x) + x - 3;
tol = fabs(p);
iterations++;
}
if(tol < 0.00001)
{
cout << "x is " << x << endl
<< iterations << " iterations ";
  
double weins_const = (h*c)/(x*k);
cout << "computed wein's_constant is: " << weins_const << endl;
cout << "real value of wein's constant is : 0.0028977729" << endl;
}
else
cout << "Did not converge after 100 iterations ";
return 0;
}

//////// Code ends here ///////////

Sample output:

x is 2.82144
3 iterations
computed wein's_constant is: 0.00509952
real value of wein's constant is : 0.0028977729

Related Questions

Navigate

• Browse (All)
• Browse W
• Subjects

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.