Gauss-eliminasjon og delvis pivotering#
TMA4400 Matematikk 1: Kalkulus og lineær algebra
Dato: 25. august 2025
Standard Gauss-eliminasjon#
La \(A\) være ei \(m\times n\)-matrise, \(\mathbf{x}\) en vektor i \(\mathbb{R}^n\) og \(\mathbf{b}\) en vektor i \(\mathbb{R}^m\). Vi skal løse ligningssystemet
ved hjelp av Gauss-eliminasjon.
Vi går ut i fra at pivot-elementene er forskjellig fra null.
import numpy as np
# Gauss-eliminasjon (antar at pivotelementene er forskjellig fra null i alle rader!)
def GE(B):
A = np.copy(B) # vi lager en kopi av den gitte matrisa B, for å unngå å endre denne
m,n = np.shape(B) # m=rader, n=kolonner; det er mulig å oppgi den utvidede matrisa, eller tom matrise
for k in range(0,m-1):
for j in range(k+1,m):
l = A[j,k]/A[k,k]
for p in range(k,n):
A[j,p] = A[j,p]-l*A[k,p]
return A
Eksempel#
La oss bruke standard Gauss-eliminasjon på følgende utvidede matrise:
A0 = np.array([[1.,-1.,1.,1.],[2.,4.,2.,1.],[1.,3.,2.,1.]]) # vi bruker 1. og så videre for å sikre flyttalltolkning i python
A1 = GE(A0)
print('\nVi begynner med \n{0}\nog bruker Gauss-eliminasjon som gitt over: \n{1}'.format(np.matrix(A0),np.matrix(A1)))
Vi begynner med
[[ 1. -1. 1. 1.]
[ 2. 4. 2. 1.]
[ 1. 3. 2. 1.]]
og bruker Gauss-eliminasjon som gitt over:
[[ 1. -1. 1. 1. ]
[ 0. 6. 0. -1. ]
[ 0. 0. 1. 0.66666667]]
Numerisk ustabilitet#
La oss bruke standard Gauss-eliminasjon på følgende utvidede matrise:
B0 = np.array([[1.e-12,1.,1.],[1.,1.,2.]])
B1 = GE(B0)
print('\nVi begynner med \n{0}\nog bruker Gauss-eliminasjon som gitt over: \n{1}'.format(np.matrix(B0),np.matrix(B1)))
Vi begynner med
[[1.e-12 1.e+00 1.e+00]
[1.e+00 1.e+00 2.e+00]]
og bruker Gauss-eliminasjon som gitt over:
[[ 1.e-12 1.e+00 1.e+00]
[ 0.e+00 -1.e+12 -1.e+12]]
Ved å regne ut for hånd får man
Her trenger vi altså en forbedret algoritme for å fange opp tilfellet hvor pivotelementet er lite, men likevel ulikt 0.
Gauss-eliminasjon med delvis pivotering#
En standard måte er såkalt delvis pivotering (engelsk: partial pivoting).
Gauss-eliminasjon med delvis pivotering skiller seg fra standard Gauss-eliminasjon ved at pivotelementet alltid er størst i absoluttverdi sammenlignet med de andre elementene i den aktuelle kolonnen under pivotelementet.
I tilfellet vi så på, altså
vil algoritmen begynne med å bytte om på rad 1 og 2 slik at vi ender opp med
Ved å gange første rad med \(-10^{-12}\) (i absoluttverdi mindre enn \(1\)) og legge denne raden til andre rad får vi:
For å illustrere algoritmen ytterligere, la oss også se på
Siden \(2>1\), må vi bytte om på rad 1 og 2 slik at vi ender opp med
Siden \(|-3|>1\), trenger vi ikke å bytte om på nåværende rad 2 og 3. Det har heller ikke noe å si at \(|-3|<4\). Vi fullfører så eliminasjonen:
# Gauss-eliminasjon med delvis pivotering
def GEPP(B):
A = np.copy(B) # vi lager en kopi av den gitte matrisa B, for å unngå å endre denne
m,n = np.shape(A) # m=rader, n=kolonner; det er mulig å oppgi den utvidede matrisa, eller tom matrise
for k in range(0,m-1):
a = np.abs(A[k,k])
pk = k
for j in range(k+1,m):
if a < np.abs(A[j,k]):
a = np.abs(A[j,k])
pk = j
if pk != k:
for p in range(k,n): # vi må bytte radene i A
r = A[k,p]
A[k,p] = A[pk,p]
A[pk,p] = r
for j in range(k+1,m):
l = A[j,k]/A[k,k]
for p in range(k,n):
A[j,p] =A[j,p]-l*A[k,p]
return A
Eksempel#
La oss bruke Gauss-eliminasjon med delvis pivotering på følgende utvidede matrise:
B2 = np.array([[1.e-12,1.,1.],[1.,1.,2.]])
B3 = GEPP(B2)
print('\nVi begynner med \n{0}\nog bruker Gauss-eliminasjon med delvis pivotering som gitt over: \n{1}'.format(np.matrix(B2),np.matrix(B3)))
Vi begynner med
[[1.e-12 1.e+00 1.e+00]
[1.e+00 1.e+00 2.e+00]]
og bruker Gauss-eliminasjon med delvis pivotering som gitt over:
[[1. 1. 2.]
[0. 1. 1.]]
Kan flere ting gå galt?#
La nå
B4 = np.array([[-1.,1.,1.,6],[1.,-1.,1.,2.],[1.,-1.,-1.,0.]])
B5 = GEPP(B4)
print('\nVi begynner med \n{0}\nog bruker Gauss-eliminasjon med delvis pivotering som gitt over: \n{1}'.format(np.matrix(B4),np.matrix(B5)))
Vi begynner med
[[-1. 1. 1. 6.]
[ 1. -1. 1. 2.]
[ 1. -1. -1. 0.]]
og bruker Gauss-eliminasjon med delvis pivotering som gitt over:
[[-1. 1. 1. 6.]
[ 0. 0. 2. 8.]
[ 0. nan nan nan]]
/tmp/ipykernel_2397/1362323269.py:18: RuntimeWarning: invalid value encountered in scalar divide
l = A[j,k]/A[k,k]
Vi deler altså på \(0\) siden pivotelementet ikke er forskjellig fra \(0\) etter hvert. Vi testet tross alt ikke for dette. La oss prøve Python sin innebygde løser:
B6 = np.array([[-1.,1.,1.],[1.,-1.,1.],[1.,-1.,-1.]])
b = np.array([6., 2., 0.])
#np.linalg.solve(B6,b) # fjern skigard-symbolet for å kjøre denne delen
La oss nå teste følgende utvidede matrise:
B7 = np.array([[1.,2.,1.,6],[1.,-1.,2.,2.],[0.8,1.,1.,0.]])
B8 = GEPP(B7)
print('\nVi begynner med \n{0}\nog bruker Gauss-eliminasjon med delvis pivotering som gitt over: \n{1}'.format(np.matrix(B7),np.matrix(B8)))
Vi begynner med
[[ 1. 2. 1. 6. ]
[ 1. -1. 2. 2. ]
[ 0.8 1. 1. 0. ]]
og bruker Gauss-eliminasjon med delvis pivotering som gitt over:
[[ 1.00000000e+00 2.00000000e+00 1.00000000e+00 6.00000000e+00]
[ 0.00000000e+00 -3.00000000e+00 1.00000000e+00 -4.00000000e+00]
[ 0.00000000e+00 0.00000000e+00 -8.32667268e-17 -4.00000000e+00]]
Utregning for hånd gir
Systemet har altså ingen løsning selv om python tror det:
B9 = np.array([[1.,2.,1.],[1.,-1.,2.],[0.8,1.,1.]])
b = np.array([6., 2., 0.])
#np.linalg.solve(B7,b) # fjern skigard-symbolet for å kjøre denne delen