Newton-Raphson
From LiteratePrograms
Contents |
The Newton-Raphson Method
This is one of the most simple methods of finding the roots of a polynomial of arbitrary degree. The method is quite easy to perform if the coefficients of the polynomial are all real and the roots of the polynomial are also real. However if the coefficients contain imaginary components, then at least one of the roots is complex and unless complex arithmetic is used, there is no way to get the complex root. After the first root is found, the coefficients of the polynomial of degree n are reduced by 1 by dividing the polynomial by the first root found. The method is then repeated until all the roots are found. This method does not involve the computation of many variables as in the case of Laguerre Method (which is based on similar principles) which it counteracts through a faster root convergence.
Quick Summary of Newton-Raphson Method
We first assume a complex root such as − 0.5 + j1. If the polynomial is given by the formula f(x), we find f( − 0.5 + j1) by direct substitution. Note that if the assumed root is the actual root, then the result of f( − 0.5 + j1) is zero because when you divide a number by its factor you get a remainder zero. After that, the polynomial is differentiated and again the assumed root is substituted to get some value. Both calculations might not yield a complex result which is not very common. The value of the f(x) is then divided by the differential value to get another complex result. This result is then subtracted from the assumed root and the process is repeated. After the root values become constant, the iteration is stopped and the polynomial is deflated with the newly found root.
Newton-Raphson Method in C#
In a new project in visual C#, add a new class. The class page will be mostly blank with the exception of some 'using directives' and the current namespace which can be changed. The best way (in my opinion) is to allow the user to input the coefficients as easily as possible. Suppose the equation : x4 + 2 * x3 + 3 * x2 + 4 * x + 5 = 0, is to be solve. The most important part of the equation is the coefficients. Therefore a most convinient way to allow the user to input such an equation is by just typing in a textbox something like: 1,2,3,4,5. Therefore all that is needed is for us to separate the coefficients and adding them to a list. In the class added, add this piece of code.
class Extract
{
public static List<double> ExtractCoeff(string input)
{
double coeff = 0;
int n = 0;
List<double> Coeff = new List<double>();
do
{
n = input.IndexOf(",");
coeff = Convert.ToDouble(input.Substring(0, n));
Coeff.Add(coeff);
input = input.Substring(n + 1);
} while (!(input.IndexOf(",") == -1));
Coeff.Add(double.Parse(input));
return (Coeff);
}
} // End ExtractCoeff Class
The code above will copy all the coefficient into a list.
Now we have to create a complex number class (assumig .net framework 4.0 is not used) which will perform operations such as subtraction, multiplication, division and exponentiation. Paste this block of code in the same class as before.
public static class Complex
{
private static List<Double> raw = new List<Double>();
public static List<double> Subtract(double x, double jx, double y, double jy)
{
raw.Clear();
raw.Add(x - y);
raw.Add(jx - jy);
return (raw);
}
public static List<double> Mult(double x, double jx, double y, double jy)
{
raw.Clear();
raw.Add(x * y - jx * jy);
raw.Add(x * jy + jx * y);
return (raw);
}
public static List<double> Div(double x, double jx, double y, double jy)
{
if (y != 0 && jy != 0)
{
raw.Clear();
double denominator = Math.Pow(y, 2) + Math.Pow(jy, 2);
Mult(x, jx, y, -jy);
raw.Add(raw[0] / denominator);
raw.Add(raw[1] / denominator);
raw.RemoveRange(0, 2);
return (raw);
}
else if (jy == 0)
{
raw.Clear();
raw.Add(x / y);
raw.Add(jx / y);
return (raw);
}
else
{
raw.Clear();
raw.Add((jx * jy) / Math.Pow(jy, 2));
raw.Add((-1 * x * jy) / Math.Pow(jy, 2));
return (raw);
}
}
public static List<double> Power(double x, double jx, double power)
{
raw.Clear();
if (jx != 0)
{
double modulus = Math.Sqrt(Math.Pow(x, 2) + Math.Pow(jx, 2));
double angle = 0;
if (x < 0 && jx > 0)
{ angle = Math.PI - Math.Atan(jx / Math.Abs(x)); }
else if (x > 0 && jx < 0)
{ angle = -Math.Atan(Math.Abs(jx) / x); }
else if (x < 0 && jx < 0)
{ angle = (Math.PI) + Math.Atan((jx) / x); }
else
{ angle = Math.Atan(jx / x); }
raw.Add(Math.Pow(modulus, power) * Math.Cos(power * angle));
raw.Add(Math.Pow(modulus, power) * Math.Sin(power * angle));
}
else
{
raw.Add(Math.Pow(x, power));
raw.Add(0);
}
return (raw);
}
} // End Complex Class
Now we're done with the complex class, we now need to find f(x) and f'(x) (first derivative). This can simply be written. We now add just below the complex class the following block of code. This class will be able to substitute a complex number into a function and differentiate and substitute. All that is needed is for the derivative arguement to be set to 1 for differentiating once.
public static class SubstDiff
{
public static List<double> SubDiff(List<double> coeff, List<double> jcoeff, double real, double imaginary, int derivative)
{
List<double> raw = new List<double>();
int degree = coeff.Count - 1; int n = 0; int s = 0;
double fx = 0; double jfx = 0;
double diffmult = 1;
for (n = 1; n <= coeff.Count - derivative - 1; n++)
{
if (derivative != 0) // This will find the number to multiply when differentiated
{
do
{
diffmult = (degree - s) * diffmult;
s++;
} while (!(s == derivative));
}
raw = Complex.Power(real, imaginary, (degree - derivative));
raw = Complex.Mult(coeff[n - 1], jcoeff[n - 1], raw[0], raw[1]);
fx += diffmult * raw[0];
jfx += diffmult * raw[1];
degree -= 1; // Reduce degree by 1
diffmult = 1; // Return it back to default
s = 0;
}
fx += coeff[n - 1];
jfx += jcoeff[n - 1];
raw.Clear();
raw.Add(fx);
raw.Add(jfx);
return (raw);
}
} // End Substitution and differentiation Class
Now we need the last class which will reduce the polynomial to a degree less when a root has been found so that we can find the other roots. One easy way to do this is by using Horner's deflation method which is straight foward. Lets assume we have the polymonial 1 * x2 − 3 * x + 2 = 0. If a root x = 2 has been found, then by Horner's method: 1. Take the coefficient of the highest degree term. In our case its 1. 2. add the coefficient of the next term (in our case -3). So: (previousvalue) * (root) + (termcoefficient). Therefore: 1 * 2 + ( − 3) = − 1
Now we've obtained two terms ; 1 and -1. Originally we had 3 coefficients, now we have two. The results can be written as ; 1 * x − 1 = 0 from which x = 1 which is actually the second root of the equation 1*x^2-3*x+2=0. You might be tempted to countinue and see what happens. Using the same formula as before "(previousvalue) * (root) + (termcoefficient)": We get; previous value = -1 ; root =2 ; term coefficient = 2; Full calculation will yield : − 1 * (2) + (2) = 0 ; which tells you that since "2" is a root, there is no remainder.
Here's the code
public static class Deflate
{
public static void Horner(ref List<double> coeff, ref List<double> jcoeff, double x, double jx)
{
int z = coeff.Count;
jcoeff.Add(jcoeff[0]);
coeff.Add(coeff[0]);
double a = 0;
double b = 0;
int c = coeff.Count;
for (int n = 1; n <= c - 3; n++)
{
List<double> raw = Complex.Mult(x, jx, coeff[coeff.Count - 1], jcoeff[coeff.Count - 1]);
a = coeff[n] + raw[0];
b = jcoeff[n] + raw[1];
coeff.Add(a);
jcoeff.Add(b);
}
coeff.RemoveRange(0, z);
jcoeff.RemoveRange(0, z);
}
} // End Deflate Class
Now we're done writing all the behind-the-scenes codes. All we have to do now is to call them. You can now use this class to create a console application (which is easier) or create a GUI. There is just very little difference in the code between the two. Below is the code written for a console. The steps are very easy to understand: 1. Assume a root 2. Get f(x) 3. Get f'(x) 4. Divide f(x) by f'(x) 5. Subtract 6. Iterate 7. Deflate and return values to default 8. Repeat all over again until there're two elements in the list 9. Add last elements and display
static void Main(string[] args)
{
List<double> coeff= new List<double>();
List<double> jcoeff= new List<double>();
List<double> roots= new List<double>();
double root = -0.5; double jroot = 1; int iterationNo = 0;
Console.WriteLine("{0}\n{1}\n{2}", "Newton-Raphson Method", "", "");
Console.WriteLine("Enter Real Coefficients");
coeff.AddRange(Extract.ExtractCoeff(Console.ReadLine()));
Console.WriteLine("Enter Imaginary Coefficients");
jcoeff.AddRange(Extract.ExtractCoeff(Console.ReadLine()));
while (coeff.Count>2)
{
while (++iterationNo <= 30)
{
List<double> rawData = new List<double>();
// Get f(x)
rawData.AddRange(SubstDiff.SubDiff(coeff, jcoeff, root, jroot, 0));
// Get f'(x)
rawData.AddRange(SubstDiff.SubDiff(coeff, jcoeff, root, jroot, 1));
// Divide f(x) by f'(x)
rawData = Complex.Div(rawData[0], rawData[1], rawData[2], rawData[3]);
// Subtract
root -= rawData[0];
jroot -= rawData[1];
rawData.Clear();
}
roots.Add(root);
roots.Add(jroot);
Deflate.Horner(ref coeff, ref jcoeff, root, jroot);
root = -0.5; jroot = 1; iterationNo = 0;
}
roots.Add(-coeff[1]);
roots.Add(-jcoeff[1]);
Console.WriteLine("");
for (int element=0;element<=roots.Count-1; element+=2)
{
Console.WriteLine("{0}\t{1}",+roots[element],"j " + roots[element + 1]);
}
Console.ReadKey();
}
Modifications
Run and test it. It should work very well except for the .net inaccuracy of 15 digits limit. This means that the more roots there is to find, the higher the degree of inaccuracy. However the gmp library for multi precision can be downloaded for free and used giving arbitrary precision in the expense of performance. Also there is another way the accuracy can be improved. There is no need to find another root when a complex root is found. By that I mean suppose you find an initial root of 1+j2. Then we know that complex numbers can never exist alone. They must have their counterpart; i.e, the conjugate pair. This means that once you get a root that is complex, then another root of similar terms exist. I mean if you get a root 1+2j, then 1-2j is another root. So, there is no need to iterate again to get that root. All you do is you deflate once with root 1+2j and again with 1-2j. This will greatly improve the speed and accuracy of the remaining roots. A degree 50 polynomial with all complex roots will no longer need to be done 50 times to get all the roots, instead 25 times will be enough. This enables such methods like Newton-Raphson and Laguerre to appear faster than others such as the Durand-Kerner method which finds all the roots simultaneously. Though the Durand-Kerner method can be more accurate than the former methods, but the facility to optimize those root-by-root methods makes them more popular for platforms where speed is a liability.
The entire project can be downloaded here [[1]]
| Download code |