Click here to Skip to main content
15,896,154 members
Articles / General Programming / Algorithms

Introduction to Numerical Solutions

Rate me:
Please Sign up or sign in to vote.
4.89/5 (11 votes)
24 Feb 2012CPOL11 min read 30.4K   970   36  
An introduction to numerical solver algorithms with general purpose demonstration code.
///////////////////////////////////////////////////////////////////////////////
//
//  Form1.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.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;

// Namespace
namespace SolverDemo
{
    /// <summary>
    /// Test class.
    /// </summary>
    public partial class Form1 : Form
    {
        /// <summary>
        /// Constructor.
        /// </summary>
        public Form1()
        {
            InitializeComponent();
        }

        /// <summary>
        /// The data to test.
        /// </summary>
        double[] data = new double[100];

        /// <summary>
        /// Compute a value.
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <param name="c"></param>
        /// <param name="x"></param>
        /// <returns></returns>
        protected double computeValue(double a, double b, double c, double x)
        {
            try
            {
                //return a * x * x + b * x * System.Math.Exp(x) + c * System.Math.Exp(x + 1);

                return System.Math.Sin(10 * a * b * c * x) +
                    System.Math.Cos(15 * a * b * x) +
                    b * x * System.Math.Log(10 + x * a / c) +
                    c * x * System.Math.Log(x + 10);
            }
            catch (Exception)
            {
                return 0;
            }
        }

        /// <summary>
        /// The error function we're using.
        /// This can be anything we can compute fast that is relative to the 
        /// error in some meaningful way.
        /// 
        /// It must be monotonic:
        ///    Reported error = F(actual error).
        ///    Where:
        ///       if actualError increases, Reported error increases
        ///       and if actual error decreases reported error decreases.
        ///
        /// In this case, the function computes the sum of the squares of the error.
        /// If we took the sqrt and divided by the number of data points
        /// we'd get the standard deviation of the error.
        /// </summary>
        /// <param name="parameters"></param>
        /// <returns></returns>
        double ComputeError(double[] parameters)
        {
            // Same as computing the data but use the guessed parameters.
            double a = parameters[0], b = parameters[1], c = parameters[2];
            double error = 0;
            for (int i = 0; i < data.Length; i++)
            {
                // di is the x in this equation so we recompute what the data is based on the 
                // parameters passed in.
                double di = (double)(i * 0.1);

                double v = computeValue(a, b, c, di);

                // v is the value if these parameters are correct, v - data[i] is the ith error.
                v = v - data[i];

                // Sum square of errors.
                error += (v * v);
            }
            return error;
        }

        



        /// <summary>
        /// Test with Paris Genetic algorithm.
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void button1_Click(object sender, EventArgs e)
        {
            Solver s = new Solver();

            // These are the "real" parameters we're trying to find.
            double a = 1.001, b = 2.002, c = 3.003;

            // Generate the fake data.
            for (int i = 0; i < data.Length; i++)
            {
                double di = (double)(i * 0.1);
                data[i] = computeValue(a, b, c, di);
            }

            // Set up each parameter's search criteria.
            Solver.Parameter pa = new Solver.Parameter(
                "a",     // Name.
                1,       // Center
                1.5,     // Max
                0.5,     // Min
                0.1,     // Initial step size
                0.00001);  // Final step size
            Solver.Parameter pb = new Solver.Parameter("b", 2, 3, 1, 0.1, 0.00001);
            Solver.Parameter pc = new Solver.Parameter("c", 3, 4, 2, 0.01, 0.00001);
            s.Parameters.Add(pa);
            s.Parameters.Add(pb);
            s.Parameters.Add(pc);

            // Solve, passing in:
            // ComputeError    The error function.
            // 10              Ten generations.
            // 1000            1k central most variations of the variables.
            // 1.5             How much over testing is done per generation.
            // 0.00005         After reaching this error stop trying.
            double err = s.Solve_ParisGenetic(ComputeError, 10, 1000, 1.25, 0.00005);

            // Our provided "error" just needed to be relative, so we never converted
            // to std. dev. to avoid many square roots. 
            //
            // Convert to stdev error.
            err = System.Math.Sqrt(err) / data.Length;

            // Output the solution.
            label1.Text = "Actual A, B, C = " + a + ", " + b + ", " + c + "\r\nSolved A, B, C = " +
                pa.Solution + ", " + pb.Solution + ", " + pc.Solution + "\r\nError = " + err;

            // Fire up the graphs
            solutionVisualizer1.Setup(s.Parameters, ComputeError, s.BestStep);
        }

        /// <summary>
        /// Solve with direct search.
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void button2_Click(object sender, EventArgs e)
        {
            Solver s = new Solver();

            // These are the "real" parameters we're trying to guess at.
            double a = 1.001, b = 2.002, c = 3.003;

            // Generate the fake data.
            for (int i = 0; i < data.Length; i++)
            {
                double di = (double)(i * 0.1);
                data[i] = computeValue(a, b, c, di);
            }

            // Set up each parameter's search criteria.
            Solver.Parameter pa = new Solver.Parameter(
                "a",     // Name.
                1,       // Center
                1.5,     // Max
                0.5,     // Min
                0.1,     // Initial step size
                0.0001);  // Final step size
            Solver.Parameter pb = new Solver.Parameter("b", 2, 3, 1, 0.1, 0.00001);
            Solver.Parameter pc = new Solver.Parameter("c", 3, 4, 2, 0.01, 0.00001);
            s.Parameters.Add(pa);
            s.Parameters.Add(pb);
            s.Parameters.Add(pc);

            // Solve:
            //   ComputeError    The error function.
            //   200             200 generations, each generation is only 10 calculations
            //   10              10 walks per variable before moving on.
            //   1.25            How much over testing is done per generation.
            //   0.00005         After reaching this error stop trying.
            double err = s.Solve_DirectSearch(ComputeError, 200, 10, 1.25, 0.00005);

            // Our provided "error" just needed to be relative, so we never converted
            // to std. dev. to avoid many square roots. 
            //
            // Convert to stdev error.
            err = System.Math.Sqrt(err) / data.Length;

            // Output results.
            label1.Text = "Actual A, B, C = " + a + ", " + b + ", " + c + "\r\nSolved A, B, C = " +
                pa.Solution + ", " + pb.Solution + ", " + pc.Solution + "\r\nError = " + err;

            // Update the graphs.
            solutionVisualizer1.Setup(s.Parameters, ComputeError, s.BestStep);
        }

        /// <summary>
        /// Slope.
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void button3_Click(object sender, EventArgs e)
        {
            Solver s = new Solver();

            // These are the "real" parameters we're trying to guess at.
            double a = 1.001, b = 2.002, c = 3.003;

            // Generate the fake data.
            for (int i = 0; i < data.Length; i++)
            {
                double di = (double)(i * 0.1);
                data[i] = computeValue(a, b, c, di);
            }

            // Set up each parameter's search criteria.
            Solver.Parameter pa = new Solver.Parameter(
                "a",     // Name.
                1.05,       // Center
                1.2,     // Max
                0.8,     // Min
                0.1,     // Initial step size
                0.0001);  // Final step size
            Solver.Parameter pb = new Solver.Parameter("b", 2.05, 2.1, 1.9, 0.1, 0.0001);
            Solver.Parameter pc = new Solver.Parameter("c", 3.05, 3.1, 2.9, 0.01, 0.0001);
            s.Parameters.Add(pa);
            s.Parameters.Add(pb);
            s.Parameters.Add(pc);

            // Pre-Solve
            //   ComputeError    The error function.
            //   200             200 generations, each generation is only 10 calculations
            //   10              10 walks per variable before moving on.
            //   1.25            How much over testing is done per generation.
            //   0.00005         After reaching this error stop trying.
            s.Solve_DirectSearch(ComputeError, 1, 10, 1.25, 0.00005);
            pa.CenterValue = pa.Solution;
            pb.CenterValue = pb.Solution;
            pc.CenterValue = pc.Solution;

            // Now re-solve with slope variate, slope variate assumes your close enough
            // to the solution that the errors can more or less fit a quadratic and the
            // trench of a sadle point can be "walked" to improve the solution.
            //
            // Thus it is here used as a refining solution.
            double err = s.Solve_SlopeVariate(ComputeError, 2, 5, 1.25, 0.00000000005);

            // Our provided "error" we can convert to stdev error.
            err = System.Math.Sqrt(err) / data.Length;

            // Output results.
            label1.Text = "Actual A, B, C = " + a + ", " + b + ", " + c + "\r\nSolved A, B, C = " +
                pa.Solution + ", " + pb.Solution + ", " + pc.Solution + "\r\nError = " + err;

            // Update graphs
            solutionVisualizer1.Setup(s.Parameters, ComputeError, s.BestStep);
        }

        /// <summary>
        /// Generation based paris genetic.
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void button4_Click(object sender, EventArgs e)
        {
            Solver s = new Solver();

            // These are the "real" parameters we're trying to find.
            double a = 1.001, b = 2.002, c = 3.003;

            // Generate the fake data.
            for (int i = 0; i < data.Length; i++)
            {
                double di = (double)(i * 0.1);
                data[i] = computeValue(a, b, c, di);
            }

            // Set up each parameter's search criteria.
            Solver.Parameter pa = new Solver.Parameter(
                "a",     // Name.
                1,       // Center
                1.5,     // Max
                0.5,     // Min
                0.1,     // Initial step size
                0.00001);  // Final step size
            Solver.Parameter pb = new Solver.Parameter("b", 2, 3, 1, 0.1, 0.00001);
            Solver.Parameter pc = new Solver.Parameter("c", 3, 4, 2, 0.01, 0.00001);
            s.Parameters.Add(pa);
            s.Parameters.Add(pb);
            s.Parameters.Add(pc);

            // Solve, passing in:
            // ComputeError    The error function.
            // 4               Four generations.
            // 1000            1k central most variations of the variables.
            // 1.5             How much over testing is done per generation.
            // 0.00005         After reaching this error stop trying.
            double err = s.Solve_ParisGeneticMultiple(ComputeError, 10, 1000, 1.25, 0.00005);

            // Our provided "error" just needed to be relative, so we never converted
            // to std. dev. to avoid many square roots. 
            //
            // Convert to stdev error.
            err = System.Math.Sqrt(err) / data.Length;

            // Output the solution.
            label1.Text = "Actual A, B, C = " + a + ", " + b + ", " + c + "\r\nSolved A, B, C = " +
                pa.Solution + ", " + pb.Solution + ", " + pc.Solution + "\r\nError = " + err;

            // Fire up the graphs
            //solutionVisualizer1.Setup(s.Parameters, ComputeError, s.BestStep);
        }
    }
}

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)


Written By
Technical Lead
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.

Comments and Discussions