Click here to Skip to main content
Click here to Skip to main content
Add your own
alternative version
Go to top

Introduction to Numerical Solutions

, 24 Feb 2012
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
{
    /// <summary>
    /// 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
    /// </summary>
    public class LinearLeastSquares
    {
        /// <summary>
        /// The solution, either:
        ///   y = Solution[0]*x + Solution[1];
        /// or if Solution.Length == 3: 
        ///   y = Solution[0] + Solution[1]*x + Solution[2]*x*x;
        /// </summary>
        public double[] Solution { get; protected set; }

        /// <summary>
        /// Solve for y = mx + b -> 
        /// y = Solution[0]*x + Solution[1];
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        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;
        }

        /// <summary>
        /// Solve for y = a0 + a1*x + a2*x^2 -> 
        /// y = Solution[0] + Solution[1]*x + Solution[2]*x*x;
        /// </summary>
        /// <param name="x"></param>
        /// <param name="y"></param>
        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;
        }

        /// <summary>
        /// Compute the estimated Y values for this x.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        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;
        }

        /// <summary>
        /// Compute the error of this fit.
        /// </summary>
        /// <param name="x"></param>
        /// <returns></returns>
        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;
        }
    }
}

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

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.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

HoshiKata
Software Developer (Senior) KMC Systems
United States 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.

| Advertise | Privacy | Mobile
Web03 | 2.8.140916.1 | Last Updated 24 Feb 2012
Article Copyright 2012 by HoshiKata
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid