14
Tags:
c#
Skrevet af
Bruger #4522
@ 22.10.2007
Crouts Faktorisering
Nu skal vi kigge på den sidste metode vi vil implementere til at løse et ligningssystem, og det er Crouts faktorisering.
De to andre metoder vi har kigget på har den ulempe, at vi ændrer b vektoren og A matricen samtidig, hvilket betyder at hvis vi vil løse et system med samme A matrix, men med forskellige b vektorer - en situation der sandsynligvis vil opstå ganske tit - skal vi enten starte forfra med reduktionen af vores A matrix eller på en eller anden måde gemme en historie over modifikationerne ved b vektoren.
I Crouts faktorisering manipulerer vi A matricen og b vektoren uafhængigt af hindanden. Vi kan derfor bruge samme A matrix sammen med flere forskellige b vektorer uden at skulle starte forfra. Denne fleksibilitet gør denne metode til nok den mest benyttede ved løsning af ligningssystemer.
Følgende ligning viser hvordan vi omdanner det oprindelige system til til et nyt system hvor A matricen er omdannet til produktet af en øvre og nedre triangulær matrix (henholdsvis U og L):
Ax = LUx = b
L og U er henholdsvis (i dette eksempel er de fra et 4x4 system; l_ij er elementet fra L matricen på i. række og j. kolonne - ligeså med u_ij for U matricens vedkommende. Desværre ligner bogstavet lille L (l) ettallet (1) ret meget, men det skulle være klart hvad der er hvad):
| 1 0 0 0 |
| l_10 1 0 0 |
L = | l_20 l_21 1 0 |
| l_30 l_31 l_32 1 |
| u_00 u_01 u_02 u_03 |
| 0 u_11 u_12 u_13 |
U = | 0 0 u_22 u_23 |
| 0 0 0 u_33 |
Indgangene i L og U udregnes kolonne for kolonne efter følgende to ligninger (a_ij er elementet i række i og kolonne j fra A matricen):
Vi overskriver A matricen (via vores a variabel) med de to triangulære matricer, og sparer derfor lidt plads på den konto. Den nye A matrix er derfor (igen bruger jeg for eksemplets skyld et 4x4 system):
| u_00 u_01 u_02 u_03 |
| l_10 u_11 u_12 u_13 |
A' = | l_20 l_21 u_22 u_23 |
| l_30 l_31 l_32 u_33 |
Med ovenstående A' matrix, løser vi følgende ligning:
Ly = b
y findes ved brug af substitution:
Når vi har fundet y, finder vi løsningsvektoren ved at løse:
Ux = y
Igen bruger vi substitution:
Vi har nu nok information til at kunne implementere Crouts løsningen i C# kode. Her er koden:
private void SolveByCrout()
{
int numberOfRows = a.Length;
int numberOfColumns = a[0].Length;
int m = 0;
double[] tempRow;
double tempElement;
// We will now get the L and U matrix as explained in the article.
// We will hold both of them in a as to save space.
for (int column = 0; column < numberOfColumns; ++column)
{
// The U matrix
for (int i = 0; i <= column; ++i)
{
for (int j = 0; j < i; ++j)
{
a[i][column] -= a[i][j] * a[j][column];
}
}
// The L matrix
for (int i = column + 1; i < numberOfRows; ++i)
{
for (int j = 0; j < column; ++j)
{
a[i][column] -= a[i][j] * a[j][column];
}
}
// The pivot
m = column;
for (int i = column + 1; i < numberOfRows; ++i)
{
if (Math.Abs(a[i][column]) > Math.Abs(a[m][column]))
{
m = i;
}
}
// New pivot: rearrange rows and the b vector
if (m != column)
{
// Rearrange a
tempRow = a[column];
a[column] = a[m];
a[m] = tempRow;
// Rearrange b
tempElement = b[column];
b[column] = b[m];
b[m] = tempElement;
}
// Divide L elements by the diagonal value
for (int i = column + 1; i < numberOfRows; ++i)
{
a[i][column] /= a[column][column];
}
}
// Now we have calculated L and U.
// Now we'll find the solution by doing the substutions as described
// in the article.
// First forward substitution
for (int i = 1; i < numberOfRows; ++i)
{
for (int j = 0; j < i; ++j)
{
b[i] -= a[i][j] * b[j];
}
}
// And now we finish with the backward substitution
b[numberOfRows - 1] = b[numberOfRows - 1] / a[numberOfRows - 1][numberOfRows - 1];
for (int i = numberOfRows - 2; i >= 0; --i)
{
for (int j = i + 1; j < numberOfColumns; ++j)
{
b[i] -= a[i][j] * b[j];
}
b[i] /= a[i][i];
}
}
Bemærk at Crouts metoden bruger sin egen "indbyggede" pivoting, hvorfor der ikke gøres brug af vores Pivot() ovenfor.
Test af EquationSolver
Til at teste EquationSolver benyttes to ligningssystemer:
Ligningssystem 1:
| 0 4 8 6 | x1 = 56
| 3 -2 2 1 | x2 = 9
| 7 -6 3 2 | x3 = 12
| 4 -3 2 4 | x4 = 20
Ligningssystem 2:
| 0 1 5 -2 | x1 = 9
| 7 -6 3 1 | x2 = 8
| 1 4 2 -3 | x3 = 3
| 5 -2 1 4 | x4 = 20
Bemærk venligst at element A[0,0] = 0 i begge eksempler, hvilket som vi nævnte under Pivoting-sektionen medfører at vores system bliver ustabilt under vores løsningsmetoder, hvorfor vi netop gør brug af pivoting.
Løsningen til testsystemerne er:
x1 = 1
x2 = 2
x3 = 3
x4 = 4
Følgende fil, EquationSolverTest.cs, indeholder unit tests for ovenstående to eksempler:
using System;
using System.Collections.Generic;
using System.Text;
using NUnit.Framework;
using NUnit.Framework.SyntaxHelpers;
namespace EquationSolver
{
[TestFixture]
public class EquationSolverTester
{
[SetUp]
public void SetUp()
{
a = new double[4][];
}
[Test]
public void TestGaussJordan()
{
a[0] = new double[] { 0, 4, 8, 6 };
a[1] = new double[] { 3, -2, 2, 1 };
a[2] = new double[] { 7, -6, 3, 2 };
a[3] = new double[] { 4, -3, 2, 4 };
double[] b = new double[] { 56, 9, 12, 20 };
EquationSolver solver = new EquationSolver(a, b, SolveTechnique.GaussJordan);
solver.Solve();
Assert.That(b[0], Is.EqualTo(1).Within(0.00000000000005));
Assert.That(b[1], Is.EqualTo(2).Within(0.00000000000005));
Assert.That(b[2], Is.EqualTo(3).Within(0.00000000000005));
Assert.That(b[3], Is.EqualTo(4).Within(0.00000000000005));
a[0] = new double[] { 0, 1, 5, -2 };
a[1] = new double[] { 7, -6, 3, 1 };
a[2] = new double[] { 1, 4, 2, -3 };
a[3] = new double[] { 5, -2, 1, 4 };
b = new double[] { 9, 8, 3, 20 };
solver = new EquationSolver(a, b, SolveTechnique.GaussJordan);
solver.Solve();
Assert.That(b[0], Is.EqualTo(1).Within(0.000000000005));
Assert.That(b[1], Is.EqualTo(2).Within(0.000000000005));
Assert.That(b[2], Is.EqualTo(3).Within(0.000000000005));
Assert.That(b[3], Is.EqualTo(4).Within(0.000000000005));
}
[Test]
public void TestGauss()
{
a[0] = new double[] { 0, 4, 8, 6 };
a[1] = new double[] { 3, -2, 2, 1 };
a[2] = new double[] { 7, -6, 3, 2 };
a[3] = new double[] { 4, -3, 2, 4 };
double[] b = new double[] { 56, 9, 12, 20 };
EquationSolver solver = new EquationSolver(a, b, SolveTechnique.Gauss);
solver.Solve();
Assert.That(b[0], Is.EqualTo(1).Within(0.00000000000005));
Assert.That(b[1], Is.EqualTo(2).Within(0.00000000000005));
Assert.That(b[2], Is.EqualTo(3).Within(0.00000000000005));
Assert.That(b[3], Is.EqualTo(4).Within(0.00000000000005));
a[0] = new double[] { 0, 1, 5, -2 };
a[1] = new double[] { 7, -6, 3, 1 };
a[2] = new double[] { 1, 4, 2, -3 };
a[3] = new double[] { 5, -2, 1, 4 };
b = new double[] { 9, 8, 3, 20 };
solver = new EquationSolver(a, b, SolveTechnique.Gauss);
solver.Solve();
Assert.That(b[0], Is.EqualTo(1).Within(0.000000000005));
Assert.That(b[1], Is.EqualTo(2).Within(0.000000000005));
Assert.That(b[2], Is.EqualTo(3).Within(0.000000000005));
Assert.That(b[3], Is.EqualTo(4).Within(0.000000000005));
}
[Test]
public void TestCrout()
{
a[0] = new double[] { 0, 4, 8, 6 };
a[1] = new double[] { 3, -2, 2, 1 };
a[2] = new double[] { 7, -6, 3, 2 };
a[3] = new double[] { 4, -3, 2, 4 };
double[] b = new double[] { 56, 9, 12, 20 };
EquationSolver solver = new EquationSolver(a, b, SolveTechnique.Crouts);
solver.Solve();
Assert.That(b[0], Is.EqualTo(1).Within(0.00000000000005));
Assert.That(b[1], Is.EqualTo(2).Within(0.00000000000005));
Assert.That(b[2], Is.EqualTo(3).Within(0.00000000000005));
Assert.That(b[3], Is.EqualTo(4).Within(0.00000000000005));
a[0] = new double[] { 0, 1, 5, -2 };
a[1] = new double[] { 7, -6, 3, 1 };
a[2] = new double[] { 1, 4, 2, -3 };
a[3] = new double[] { 5, -2, 1, 4 };
b = new double[] { 9, 8, 3, 20 };
solver = new EquationSolver(a, b, SolveTechnique.Crouts);
solver.Solve();
Assert.That(b[0], Is.EqualTo(1).Within(0.000000000005));
Assert.That(b[1], Is.EqualTo(2).Within(0.000000000005));
Assert.That(b[2], Is.EqualTo(3).Within(0.000000000005));
Assert.That(b[3], Is.EqualTo(4).Within(0.000000000005));
}
double[][] a;
}
}
Ovenstående unit tests viser at vores ligningsløser fungerer rigtig godt! For at være helt sikker skal man selvfølgelig kører implementationen igennem mange flere tests for at sikre at koden holder i alle tilfælde, inkl. grænsetilfælde osv.
Bemærk at vi arbejder med double's hvilket vil sige at vi skal bruge en tolerance parameter ved testen som angiver hvor meget resultatværdien må afvige fra den forventede værdi. Jeg har anvendt en meget lille toleranceværdi (0.000000000005) som tit vil være påkrævet ved videnskabelige udregninger. Til finans/business og andre formål kan en toleranceværdi på f.eks. 0.0005 tit vær god nok.
Så! Nu er vi færdig. I denne artikel kiggede vi på hvordan vi i C#/.NET kunne implementere de mest almindelige løsningsmetoder for løsning af ligningssystemer. Jeg håber det har været ligeså interessant at læse som det har været at skrive.
Hvad synes du om denne artikel? Giv din mening til kende ved at stemme via pilene til venstre og/eller lægge en kommentar herunder.
Del også gerne artiklen med dine Facebook venner:
Kommentarer (3)
Yderst interessant artikel! Flere af den slags
5 herfra!
Det er fedt med nogle mere avancerede artikler her på Udvikleren, jeg kunne dog ikke helt følge med i Crouts Faktorisering, måske pga. det sene tidspunkt.
Det er lidt uheldigt at det eksempel på et ligningssystem du angiver:
x² + y² = 1
2x+ 4y = 0
Ikke er lineært (pga. x² og y²) og derved ikke kan repræsenteres af en matrix og altså ikke løses vha. disse fremgangsmåder.
Du burde også læse om C#'s rektangulære arrays, en af C#'s fordele i forhold til Java. Der er en rigtig god illustration her
http://decenturl.com/books.google/rectangular-vs-jagged og her er der en god gennemgang af fordele ulemper ved de 2 måder
http://msdn.microsoft.com/msdnmag/issues/04/03/ScientificC/#S10
Hej Martin,
Tak for dine kommentarer.
Jeg bemærkede slet ikke at mit eksempel på et ligningssytem ikke var lineært, hvilket ikke er særlig smart. Tak fordi du gjorde mig opmærksom på det - det vil jeg få ændret.
Og tak for linksene på C#s rektangulære arrays, dem vil jeg læse.
Du skal være
logget ind for at skrive en kommentar.