locked
Solution of linear equations systems

    Question

  • Hi,
    I need a simple function that will solve a system of linear equations. My problems are in the order of 10x10 equations, so I need something very simple. Any idea of where I could find one of these functions?

    Thanks for your help,
    Mike
    Friday, January 22, 2010 7:28 PM

Answers

  • As with most numerical routines, it is best to use a library and not code this yourself as there are pitfalls in ensuring the accuracy of the computations.

    There is a list of libraries, some open source, here: http://en.wikipedia.org/wiki/List_of_numerical_libraries
    • Marked as answer by Harry Zhu Friday, January 29, 2010 2:42 AM
    Saturday, January 23, 2010 12:00 AM
  • Use Guassian Elimination.  I've got source code on these forums somwhere.

    This function takes an N X N matrix in the X variable and a 1 X N vector in the Y variable.  The coefficients are returned in the Y variable.

        public void ComputeCoefficents(double[,] X, double[] Y)
        {
          
    int I, J, K, K1, N;
          N = Y.Length;
          
    for (K = 0; K < N; K++)
          {
            K1 = K + 1;
            
    for (I = K; I < N; I++)
            {
              
    if (X[I, K] != 0)
              {
                
    for (J = K1; J < N; J++)
                {
                  X[I, J] /= X[I, K];
                }
                Y[I] /= X[I, K];
              }
            }
            
    for (I = K1; I < N; I++)
            {
              
    if (X[I, K] != 0)
              {
                
    for (J = K1; J < N; J++)
                {
                  X[I, J] -= X[K, J];
                }
                Y[I] -= Y[K];
              }
            }
          }
          
    for (I = N - 2; I >= 0; I--)
          {
            
    for (J = N - 1; J >= I + 1; J--)
            {
              Y[I] -= X[I, J] * Y[J];
            }
          }
        }
      }
    • Marked as answer by Harry Zhu Friday, January 29, 2010 2:42 AM
    Friday, January 22, 2010 11:53 PM
  • Just in case, here is a function I had to write myself. It includes a sample of it usage and a check of the results.

            static void Main(string[] args)
            {
                double[,] cc = {{2,1,-1},{-3,-1,2},{-2,1,2}};
                double[,] dd = new double[3, 3];
                for (int ii=0;ii<3;ii++)
                {
                    for (int jj = 0; jj < 3; jj++)
                    {
                        dd[ii,jj] = cc[ii,jj];
                    }
                }
                double[]xxx={8,-11,-3};
                double[]rrx=new double[3];
                int nn = 3;
                linEq(cc, xxx, rrx, nn);
                double[] rry = new double[3];
                for (int ii = 0; ii < 3; ii++)
                {
                    rry[ii] = 0;
                    for (int jj = 0; jj < 3; jj++)
                    {
                        rry[ii] += dd[ii, jj] * rrx[jj];
                    }
                }
            }
            static void linEq(double[,]cc,double[]xx,double[] rr,int nn)
            {
                // forward elimination
                for (int ii = 0; ii < nn-1; ii++)
                {
                    // pivot row
                    double pp = cc[ii,ii];
                    for (int kk = ii; kk < nn; kk++)
                    {
                        cc[ii, kk] /= pp;
                    }
                    xx[ii] /= pp;
                    // eliminate unknown ii from rows below
                    for (int jj = ii+1; jj < nn; jj++)
                    {
                        double dd = cc[jj, ii];
                        for (int kk = ii; kk < nn; kk++)
                        {
                            cc[jj, kk] -= cc[ii, kk] * dd;
                        }
                        xx[jj] -= xx[ii] * dd;
                    }
                }

                // backward substitution
                rr[nn - 1] = xx[nn - 1] / cc[nn - 1, nn - 1];
                for (int ii = nn - 2; ii >= 0; ii--)
                {
                    rr[ii] = xx[ii];
                    for (int jj = ii+1; jj < nn; jj++)
                    {
                        rr[ii] -= cc[ii,jj] * rr[jj];
                    }
                }
            }

    • Marked as answer by Harry Zhu Friday, January 29, 2010 2:42 AM
    Tuesday, January 26, 2010 10:45 PM

All replies

  • Use Guassian Elimination.  I've got source code on these forums somwhere.

    This function takes an N X N matrix in the X variable and a 1 X N vector in the Y variable.  The coefficients are returned in the Y variable.

        public void ComputeCoefficents(double[,] X, double[] Y)
        {
          
    int I, J, K, K1, N;
          N = Y.Length;
          
    for (K = 0; K < N; K++)
          {
            K1 = K + 1;
            
    for (I = K; I < N; I++)
            {
              
    if (X[I, K] != 0)
              {
                
    for (J = K1; J < N; J++)
                {
                  X[I, J] /= X[I, K];
                }
                Y[I] /= X[I, K];
              }
            }
            
    for (I = K1; I < N; I++)
            {
              
    if (X[I, K] != 0)
              {
                
    for (J = K1; J < N; J++)
                {
                  X[I, J] -= X[K, J];
                }
                Y[I] -= Y[K];
              }
            }
          }
          
    for (I = N - 2; I >= 0; I--)
          {
            
    for (J = N - 1; J >= I + 1; J--)
            {
              Y[I] -= X[I, J] * Y[J];
            }
          }
        }
      }
    • Marked as answer by Harry Zhu Friday, January 29, 2010 2:42 AM
    Friday, January 22, 2010 11:53 PM
  • As with most numerical routines, it is best to use a library and not code this yourself as there are pitfalls in ensuring the accuracy of the computations.

    There is a list of libraries, some open source, here: http://en.wikipedia.org/wiki/List_of_numerical_libraries
    • Marked as answer by Harry Zhu Friday, January 29, 2010 2:42 AM
    Saturday, January 23, 2010 12:00 AM
  • Just in case, here is a function I had to write myself. It includes a sample of it usage and a check of the results.

            static void Main(string[] args)
            {
                double[,] cc = {{2,1,-1},{-3,-1,2},{-2,1,2}};
                double[,] dd = new double[3, 3];
                for (int ii=0;ii<3;ii++)
                {
                    for (int jj = 0; jj < 3; jj++)
                    {
                        dd[ii,jj] = cc[ii,jj];
                    }
                }
                double[]xxx={8,-11,-3};
                double[]rrx=new double[3];
                int nn = 3;
                linEq(cc, xxx, rrx, nn);
                double[] rry = new double[3];
                for (int ii = 0; ii < 3; ii++)
                {
                    rry[ii] = 0;
                    for (int jj = 0; jj < 3; jj++)
                    {
                        rry[ii] += dd[ii, jj] * rrx[jj];
                    }
                }
            }
            static void linEq(double[,]cc,double[]xx,double[] rr,int nn)
            {
                // forward elimination
                for (int ii = 0; ii < nn-1; ii++)
                {
                    // pivot row
                    double pp = cc[ii,ii];
                    for (int kk = ii; kk < nn; kk++)
                    {
                        cc[ii, kk] /= pp;
                    }
                    xx[ii] /= pp;
                    // eliminate unknown ii from rows below
                    for (int jj = ii+1; jj < nn; jj++)
                    {
                        double dd = cc[jj, ii];
                        for (int kk = ii; kk < nn; kk++)
                        {
                            cc[jj, kk] -= cc[ii, kk] * dd;
                        }
                        xx[jj] -= xx[ii] * dd;
                    }
                }

                // backward substitution
                rr[nn - 1] = xx[nn - 1] / cc[nn - 1, nn - 1];
                for (int ii = nn - 2; ii >= 0; ii--)
                {
                    rr[ii] = xx[ii];
                    for (int jj = ii+1; jj < nn; jj++)
                    {
                        rr[ii] -= cc[ii,jj] * rr[jj];
                    }
                }
            }

    • Marked as answer by Harry Zhu Friday, January 29, 2010 2:42 AM
    Tuesday, January 26, 2010 10:45 PM