12,888,335 members (45,716 online)

#### Stats

17.2K views
34 bookmarked
Posted 24 Feb 2012

# Introduction to Numerical Solutions

, 24 Feb 2012 CPOL
An introduction to numerical solver algorithms with general purpose demonstration code.
 ```﻿/////////////////////////////////////////////////////////////////////////////// // // Program.cs // // By Philip R. Braica (HoshiKata@aol.com, VeryMadSci@gmail.com) // // Distributed under the The Code Project Open License (CPOL) // http://www.codeproject.com/info/cpol10.aspx /////////////////////////////////////////////////////////////////////////////// // Using. using System; using System.Collections.Generic; using System.Linq; using System.Text; // Namespace. namespace SolverDemo { /// /// Linear least squares for 2 variables, 1st and second order, /// beyond that requires a decent matrix math package with good pseudo-inverse. /// Note that non-linear is more accurate when there is noise but overkill for /// this application, this is just a tiny quick fitter for use by the Solver class. /// /// See the following for more information: /// http://en.wikipedia.org/wiki/Simple_linear_regression /// http://en.wikipedia.org/wiki/Invertible_matrix /// public class LinearLeastSquares { /// /// The solution, either: /// y = Solution[0]*x + Solution[1]; /// or if Solution.Length == 3: /// y = Solution[0] + Solution[1]*x + Solution[2]*x*x; /// public double[] Solution { get; protected set; } /// /// Solve for y = mx + b -> /// y = Solution[0]*x + Solution[1]; /// /// /// public bool Solve_FirstOrder(double[] x, double[] y) { // This is the classical software approach to doing // this fast, form the sums, do as little math as needed. double sx = 0; double sxx = 0; double sy = 0; double sxy = 0; double syy = 0; // compute the length. int n = x.Length < y.Length ? x.Length : y.Length; // Sum. for (int i = 0; i < n; i++) { double xi = x[i]; double yi = y[i]; sx += xi; sxx += xi * xi; sy += yi; sxy += xi * yi; syy += yi * yi; } // y = bx + a double b = ((n * sxy) - (sx * sy)) / ((n * sxx) - (sx * sx)); double a = (sy - (b*sx))/n; // Store. Solution = new double[2]; Solution[0] = a; Solution[1] = b; // Done. return true; } /// /// Solve for y = a0 + a1*x + a2*x^2 -> /// y = Solution[0] + Solution[1]*x + Solution[2]*x*x; /// /// /// public bool Solve_SecondOrder(double[] x, double[] y) { // This is not done as a matrix operation // because it is much faster to calculate this way in part // because we know of the inherent matrix symetry and 3x3 // matrixes aren't soo bad to invert in code this way. // // Each is a sum: // sx = sum of the x values. // sx2 = sum of each x value after squaring that value. double sx = 0; double sx2 = 0; double sx3 = 0; double sx4 = 0; double sy = 0; double sxy = 0; double sx2y = 0; double sy2 = 0; // n is the data set length. int n = x.Length < y.Length ? x.Length : y.Length; // Sum everything. for (int i = 0; i < n; i++) { // Avoid repeated de-indexing, and repeated multiplies. double xi = x[i]; double yi = y[i]; double xi2 = xi * xi; // Sum. sx += xi; sx2 += xi2; sy += yi; sxy += xi * yi; sy2 += yi * yi; sx3 += xi2 * xi; sx4 += xi2 * xi2; sx2y += xi2 * yi; } // In matrix form, solving for [a0, a1, a2]. // n, sx, sx2 a0 sy // sx, sx2, sx3 * a1 = sxy // sx2, sx3, sx4 a2 sx2y // // Has a solution if the determinant isn't zero of: // a, b, c sy a0 // b, d, e * sxy * 1/det = a1 // c, e, f sx2y a2 // det = n * a + sx * b + sx2 * c double a = (sx2 * sx4) - (sx3 * sx3); double b = (sx3 * sx2) - (sx * sx4); double c = (sx * sx3) - (sx2 * sx2); double d = (n * sx4) - (sx2 * sx2); double e = (sx * sx2) - (n * sx3); double f = (n * sx2) - (sx * sx); double det = (n * a) + (sx * b) + (sx2 * c); // Don't bother trying if (det == 0) return false; det = 1 / det; double a0 = det * ((sy * a) + (sxy * b) + (sx2y * c)); double a1 = det * ((sy * b) + (sxy * d) + (sx2y * e)); double a2 = det * ((sy * c) + (sxy * e) + (sx2y * f)); // Store Solution = new double[3]; Solution[0] = a0; Solution[1] = a1; Solution[2] = a2; // Done. return true; } /// /// Compute the estimated Y values for this x. /// /// /// public double[] ComputeYValues(double[] x) { double[] y = new double[x.Length]; if (Solution.Length == 2) { double b = Solution[0]; double a = Solution[1]; for (int i = 0; i < x.Length; i++) { y[i] = (a * x[i]) + b; } return y; } if (Solution.Length == 3) { double a0 = Solution[0]; double a1 = Solution[1]; double a2 = Solution[2]; for (int i = 0; i < x.Length; i++) { double xi = x[i]; y[i] = a0 + (a1 * xi) + (a2 * xi * xi); } return y; } return null; } /// /// Compute the error of this fit. /// /// /// public double ComputeError(double[] x, double[] y) { double[] ye = ComputeYValues(x); double err = 0; for (int i = 0; i < x.Length; i++) { double de = y[i] - ye[i]; err += (de * de); } err = System.Math.Sqrt(err) / x.Length; return err; } } } ```

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

## Share

 Software Developer (Senior) KMC Systems United States
Phil is a Principal Software developer focusing on weird yet practical algorithms that run the gamut of embedded and desktop (PID loops, Kalman filters, FFTs, client-server SOAP bindings, ASIC design, communication protocols, game engines, robotics).

In his personal life he is a part time mad scientist, full time dad, and studies small circle jujitsu, plays guitar and piano.

 Pro