Solve a system of equations with Gaussian elimination in C#

example

This example shows how to use Gaussian elimination to solve a linear system of equations of the form:

    A1*x1 + B1*x2 + ... + N1*xn = C1
    A2*x1 + B2*x2 + ... + N2*xn = C2
    ...
    An*x1 + Bn*x2 + ... + Nn*xn = Cn

For example:

    9 * x1 + 4 * x2 = 7
    4 * x1 + 3 * x2 = 8

The solution to these equations gives values for x1 and x2 that make both equations true. (In this example, you can verify that the values -1 and 4 work.)

You can represent these equations as a matrix multiplied by a vector of variables x1, x2, …, xn, resulting in a vector of constants C1, C2, …, Cn.

    | A1  B1  ...  N1 |   | x1 |   | C1 |
    | A2  B2  ...  N2 |   | x2 |   | C2 |
    |         ...     | * |... | = |... |
    | An  Bn  ...  Nn |   | xn |   | Cn |

The helper method LoadArray loads data from text boxes into an augmented matrix that includes a column holding the constants and a final column that will eventually hold the solution. The augmented array has the form:

    | A1  B1  ...  N1  C1 0 |
    | A2  B2  ...  N2  C2 0 |
    |         ...         0 |
    | Am  Bm  ...  Nm  Cm 0 |

The following code shows the LoadArray method.

// Load the augmented array.
// Column num_cols holds the result values.
// Column num_cols + 1 will hold the variables'
// final values after backsolving.
private double[,] LoadArray(out int num_rows, out int num_cols)
{
    // Build the augmented matrix.
    string[] value_rows = txtValues.Text.Split(
        new string[] {"\r\n"},
        StringSplitOptions.RemoveEmptyEntries);
    string[] coef_rows = txtCoefficients.Text.Split(
        new string[] { "\r\n" },
        StringSplitOptions.RemoveEmptyEntries);
    string[] one_row = coef_rows[0].Split(
        new string[] { " " },
        StringSplitOptions.RemoveEmptyEntries);
    num_rows = coef_rows.GetUpperBound(0) + 1;
    num_cols = one_row.GetUpperBound(0) + 1;
    double[,] arr = new double[num_rows, num_cols + 2];
    for (int r = 0; r < num_rows; r++)
    {
        one_row = coef_rows[r].Split(
            new char[] {' '},
            StringSplitOptions.RemoveEmptyEntries);
        for (int c = 0; c < num_cols; c++)
        {
            arr[r, c] = double.Parse(one_row[c]);
        }
        arr[r, num_cols] = double.Parse(value_rows[r]);
    }

    return arr;
}

When you click the Solve button, the program uses row operations to zero out leading terms in every row except one. For example, to zero out the leading term in the second row, the program multiplies each entry in the first row by -A2 / A1 and adds the result to the second row. Multiplying A1 by -A2 / A1 gives -A2, which cancels out the value in the first column of row 2. The program makes similar substitutions to zero out the first element in the other rows.

Next the program uses row 2 to zero out the second column in rows 3 and later. It continues in this way until it can no longer zero out rows. Ideally at that point the final row has the form | 0 0 ... Kn Ln |. This represents the equation 0*x1 + 0*x2 + ... + Kn*xn = Ln so it is easy to solve for xn. In that case, xn = Ln / Kn. The program plugs that value into the second-to-last row to calculate x(n-1) and continues “backsolving” up the list until it has a value for every variable.

// Solve the system of equations.
private void cmdSolve_Click(object sender, EventArgs e)
{
    const double tiny = 0.00001;
    string txt = "";

    // Build the augmented matrix.
    // The values num_rows and num_cols are the number of rows
    // and columns in the matrix, not the augmented matrix.
    int num_rows, num_cols;
    double[,] arr = LoadArray(out num_rows, out num_cols);
    double[,] orig_arr = LoadArray(out num_rows, out num_cols);
    
    // Display the initial arrays.
    PrintArray(arr);
    PrintArray(orig_arr);

    // Start solving.
    for (int r = 0; r < num_rows - 1; r++)
    {
        // Zero out all entries in column r after this row.
        // See if this row has a non-zero entry in column r.
        if (Math.Abs(arr[r, r]) < tiny)
        {
            // Too close to zero. Try to swap with a later row.
            for (int r2 = r + 1; r2 < num_rows; r2++)
            {
                if (Math.Abs(arr[r2, r]) > tiny)
                {
                    // This row will work. Swap them.
                    for (int c = 0; c <= num_cols; c++)
                    {
                        double tmp = arr[r, c];
                        arr[r, c] = arr[r2, c];
                        arr[r2, c] = tmp;
                    }
                    break;
                }
            }
        }

        // If this row has a non-zero entry in column r, use it.
        if (Math.Abs(arr[r, r]) > tiny)
        {
            // Zero out this column in later rows.
            for (int r2 = r + 1; r2 < num_rows; r2++)
            {
                double factor = -arr[r2, r] / arr[r, r];
                for (int c = r; c <= num_cols; c++)
                {
                    arr[r2, c] = arr[r2, c] + factor * arr[r, c];
                }
            }
        }
    }

    // Display the upper-triangular array.
    PrintArray(arr);

    // See if we have a solution.
    if (arr[num_rows - 1, num_cols - 1] == 0)
    {
        // We have no solution.
        // See if all of the entries in this row are 0.
        bool all_zeros = true;
        for (int c = 0; c <= num_cols + 1; c++)
        {
            if (arr[num_rows - 1, c] != 0)
            {
                all_zeros = false;
                break;
            }
        }
        if (all_zeros)
        {
            txt = "The solution is not unique";
        }
        else
        {
            txt = "There is no solution";
        }
    }
    else
    {
        // Backsolve.
        for (int r = num_rows - 1; r >= 0; r--)
        {
            double tmp = arr[r, num_cols];
            for (int r2 = r + 1; r2 < num_rows; r2++)
            {
                tmp -= arr[r, r2] * arr[r2, num_cols + 1];
            }
            arr[r, num_cols + 1] = tmp / arr[r, r];
        }

        // Display the results.
        txt = "       Values:";
        for (int r = 0; r < num_rows; r++)
        {
            txt += "\r\nx" + r.ToString() + " = " +
                arr[r, num_cols + 1].ToString();
        }

        // Verify.
        txt += "\r\n    Check:";
        for (int r = 0; r < num_rows; r++)
        {
            double tmp = 0;
            for (int c = 0; c < num_cols; c++)
            {
                tmp += orig_arr[r, c] * arr[c, num_cols + 1];
            }
            txt += "\r\n" + tmp.ToString();
        }

        txt = txt.Substring("\r\n".Length + 1);
    }

    txtResults.Text = txt;
}

There’s one problem with this method. Suppose you are trying to zero out items in column K but the value in position [K, K] is 0. Then you can’t divide by the value to zero out other entries in that column. In that case, you can swap this row with one of the following rows that has a non-zero entry in that column and continue.

If there is no row to swap with, there are two possibilities. First, if every remaining entry in the array is zero, then the previous rows are true no matter what values the the following variables have so there is an infinite number of solutions. For example, consider the equations 1 * x0 - 1 * x1 = 2 and 2 * x0 - 2 * x1 = 4. These equations are true if x0 = x1 no matter what value x0 and x1 have.

The second possibility if you cannot finish zeroing out columns is that the final rows contain non-zero entries. In that case there is no solution possible. For example, the equations x0 + x1 = 0 and x0 + x1 = 1 do not have a solution.

Here are some other equations you can try to test the program:

| 1 -3   1|  |x1|   | 4|
| 2 -8   8| *|x2| = |-2| has solution 3, -1, -2.
|-6  3 -15|  |x3|   | 9|


| 9  3  4|  |x1|   | 7|
| 4  3  4| *|x2| = | 8| has solution -0.2, 4, -0.8.
| 1  1  1|  |x3|   | 3|


|   1   1   1  1 1|   |x1|   |  1|
|  32  16   8  4 2|   |x2|   | -1|
| 243  81  27  9 3| * |x3| = |  8|
|1024 256  64 16 4|   |x4|   |-56|
|3125 625 125 25 5|   |x5|   |569|

Has solution:
 x0 = 7.86666666666713
 x1 = -82.7500000000049
 x2 = 302.166666666684
 x3 = -446.750000000026
 x4 = 220.466666666679


| 2 -1  1|  |x1|   |1|
| 3  2 -4| *|x2| = |4| has no solution.
|-6  3 -3|  |x3|   |2|


| 1 -1  2|  |x1|   |-3|
| 4  4 -2| *|x2| = | 1| has no unique solution.
|-2  2 -4|  |x3|   | 6|

For a more in-depth discussion of Gaussian elimination, see my article Predicting Your Firm’s Future with Least Squares, Part II.

For a more general and theoretical discussion on Gaussian elimination, see the article Gaussian Elimination by Eric W. Weisstein at MathWorld–A Wolfram Web Resource.


Download Example   Follow me on Twitter   RSS feed   Donate




This entry was posted in algorithms, mathematics and tagged , , , , , , , , , . Bookmark the permalink.

5 Responses to Solve a system of equations with Gaussian elimination in C#

  1. Pingback: Select a conic section in C# - C# HelperC# Helper

  2. Mike says:

    Did you even test this?
    LoadArray is loading the same first coefficient input row for all rows of the array..

    • RodStephens says:

      Did you? I did before and I just did again, and it seem to work. Did you modify it? Or perhaps I’m not understanding your comment. The program displays the array after calling LoadArray in the console window and lot looks correct to me.

  3. Hi everyone 🙂 I’m rather new to 3d printing and I have many questions on the topic, so I hope you will not get mad at me for asking here at least couple of them. I think before I’ll get seriously into designing I should focus on the software itself, and that’s what I would like to ask you about. Mainly, should I begin with the most simple CAD there is or would it be better to start on something more complicated? I’m worried that I’ll get some undesirable habits while working on less complex software. My second question is about the software as well: should I search for CAD software that would allow me design and slice it in it, or should I use a different program for each of them? Does it even make a difference? Surprisingly, I couldn’t find the answer to that, as it seems like most websites want to focus on the very basics (like what is 3d printing and so on), and while the answers to those questions are fine, it seems like no one wants to go into the details (it looks like some of them even steal from each other! I swear I’ve found the same answers to the same questions on at least 3 different articles) but I’m getting off-topic… The last question is about 3d pens. Would it be possible to somehow convert whatever I draw with a 3d pen to a 3d model in a program? For example, if I’ll draw a horse with 3d pen, would it be possible to get its design in a program? I’m not sure how that could even work, but the very idea sounds interesting to me. Anyway, I think I’ll stop here just in case no one will ever answer me and all of this writing will be for nothing. I’m sorry that I’m using your content to ask questions, but I hope you can relate and assist a rookie like me. Anyway, thank you for posting. I did learn something from this and that’s always appreciated. Thank you, and I hope to hear back from you very soon 🙂

    • RodStephens says:

      I don’t have a 3D printer, so I can’t help much. The only question that I think I can answer is the one about the 3D pen. I think those basically just squirt out a thread of plastic so they don’t have any memory of where they are putting it. That means they don’t have any sort of model that you could transfer into a printer. At best you would need to get a 3D scanning system to scan the horse or whatever.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.