Linear Equations in IronPython QuickStart Sample

Illustrates how to solve systems of simultaneous linear equations in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to solve systems of linear equations using Numerics.NET. It shows both basic solution methods and more advanced techniques using LU decomposition.

The sample illustrates several key concepts:

• Creating and populating matrices with data
• Solving a system of equations with a single right-hand side using the Solve method
• Solving systems with multiple right-hand sides efficiently
• Checking matrix properties like singularity and condition number
• Computing matrix inverses and determinants
• Using LU decomposition for improved performance with multiple operations
• Accessing the L and U factors of the decomposition

Through clear examples and detailed output, you’ll learn the recommended approaches for solving linear systems and understanding their numerical properties. The sample also demonstrates best practices for handling both single and multiple right-hand sides, and shows how to evaluate the quality of solutions through condition number estimation.

The code

import numerics

# The DenseMatrix and LUDecomposition classes reside in the 
# Extreme.Mathematics.LinearAlgebra namespace.
from Extreme.Mathematics import *
from Extreme.Mathematics.LinearAlgebra import *

# Illustrates solving systems of simultaneous linear
# equations using the DenseMatrix and LUDecomposition classes 
# in the Extreme.Mathematics.LinearAlgebra namespace of the Extreme 
# Optimization Mathematics Library for .NET.

# A system of simultaneous linear equations is
# defined by a square matrix A and a right-hand
# side B, which can be a vector or a matrix.
#
# You can use any matrix type for the matrix A.
# The optimal algorithm is automatically selected.

# Let's start with a general matrix:
m = Matrix([[1, 1, 1, 1], [1, 2, 4, 2], [1, 3, 9, 1], [1, 4, 16, 2]])
b1 = Vector([ 1, 3, 6, 3 ])
b2 = Matrix([[1, 6], [2, 5], [3, 3], [3, 8]])
print "m = {0:.4f}".format(m)
			
#
# The Solve method
#

# The following solves m x = b1. The second 
# parameter specifies whether to overwrite the
# right-hand side with the result.
x1 = m.Solve(b1, False)
print "x1 = {0:.4f}".format(x1)
# If the overwrite parameter is omitted, the
# right-hand-side is overwritten with the solution:
m.Solve(b1)
print "b1 = {0:.4f}".format(b1)
# You can solve for multiple right hand side 
# vectors by passing them in a DenseMatrix:
x2 = m.Solve(b2, False)
print "x2 = {0:.4f}".format(x2)
			
#
# Related Methods
#

# You can verify whether a matrix is singular
# using the IsSingular method:
print "IsSingular(m) =", m.IsSingular()
# The inverse matrix is returned by the Inverse
# method:
print "Inverse(m) = {0:.4f}".format(m.GetInverse())
# The determinant is also available:
print "Det(m) = {0:.4f}".format(m.GetDeterminant())
# The condition number is an estimate for the
# loss of precision in solving the equations
print "Cond(m) = {0:.4f}".format(m.EstimateConditionNumber())
print 

#
# The LUDecomposition class
#

# If multiple operations need to be performed
# on the same matrix, it is more efficient to use
# the LUDecomposition class. This class does the
# bulk of the calculations only once.
print "Using LU Decomposition:"
# The constructor takes an optional second argument
# indicating whether to overwrite the original
# matrix with its decomposition:
lu = m.GetLUDecomposition(False)
# All methods mentioned earlier are still available:
x2 = lu.Solve(b2, False)
print "x2 = {0:.4f}".format(x2)
print "IsSingular(m) =", lu.IsSingular()
print "Inverse(m) = {0:.4f}".format(lu.GetInverse())
print "Det(m) = {0:.4f}".format(lu.GetDeterminant())
print "Cond(m) = {0:.4f}".format(lu.EstimateConditionNumber())
# In addition, you have access to the
# components, L and U of the decomposition.
# L is lower unit-triangular:
print "  L = {0:.4f}".format(lu.LowerTriangularFactor)
# U is upper triangular:
print "  U = {0:.4f}".format(lu.UpperTriangularFactor)