Øving 4

• 6.1: 1, 12, 31, 36 (scan)
• 6.2: 3, 13, 26 (scan)

Se siden for gamle eksamensoppgaver.

• Sommer 2008: 6
• Høst 2007: 1a

Oppgaven inneholder en masse tekst. Nesten alt er hjelpetekst, og selve oppgaven er svært kort. Fatt mot!

LU-faktorisering er av grunnleggende viktighet i lineæralgebra, så programvare for numerikk (som NumPy/SciPy, Matlab og LAPACK) inneholder ferdigskrevne rutiner for å beregne faktoriseringen (og for basert på den å finne løsning av ligningssystemer). Denne oppgaven vil demonstrere SciPys LU-faktorisering.

SciPys LU-faktorisering beregner LU-faktorisering med pivotering, altså en faktorisering på formen \(A=PLU\), hvor \(P\) bare er en permutasjonsmatrise. Som dokumentasjonen tilsier returnerer lu_factor ikke de tre matrisene \(P\), \(L\) og \(U\) direkte, men en mer kompakt form som forklart under.

Følgende eksempel viser hvordan vi kan faktorisere matrisen fra oppgave 20.2.4 på øving 3:

# -*- coding: utf-8 -*-
from scipy.linalg import lu_factor, lu_solve
import numpy as np
 
a = np.array([[2, 1, 2],
              [-2, 2, 1],
              [1, 2, -2]])
 
l_og_u, p = lu_factor(a)  # lu_factor returnerer to verdier.
                          # Den første, her kalt l_og_u har U-matrisen på diagonalen og
                          # i øvre triangel, og L-matrisen i nedre triangel.
                          # l_og_u er altså U pluss L minus enhetsmatrisen.
                          # p er en vektor som uttrykker permutasjonsmatrisen: Hvis element
                          # i har verdi j, vil det si at rad i og j ble byttet om i faktoriseringen
                          # av a.
 
print l_og_u
print p

Ved kjøring gir dette

[[ 2.   1.   2. ]
 [-1.   3.   3. ]
 [ 0.5  0.5 -4.5]]
[0 1 2]

som forventet (som vi ser er \(P\) her bare identitetsmatrisen).

lu_solve lar en løse ligningssystemer gitt en LU-faktorisering fra lu_solve. Vi kan utvide koden over med følgende for å løse ligningssystemet fra nevnte oppgave:

b = np.array([[0],
              [0],
              [18]])
 
print "Løsning:"
print lu_solve((l_og_u, p), b)

Dette gir:

Løsning:
[[ 2.]
 [ 4.]
 [-4.]]

1. Løs oppgave 20.2.3 ved hjelp av lu_solve og lu_factor. Altså: Gjør akkurat som over for ligningssystemet i oppgave 20.2.3.
2. Hvis du har tid: Løs 20.2.11 ved hjelp av de tilsvarende cho_factor og cho_solve for Cholesky-faktorisering.

I denne oppgaven skal ikke noe programmeres på egenhånd. Følgende Python-kode er gitt (du trenger ikke kjøre koden, men kan selvfølgelig gjøre det om du vil):

import numpy as np
 
def foo(a, b, x0, eps):
    n = x0.shape[0]       # Lengden til vektoren x0.
    x = np.zeros(n)
    for j in range(0, n):
        x[j] = x0[j]
 
    done = False
 
    while not done:
        xprev = np.zeros(n)
        for j in range(0, n):
            xprev[j] = x[j]
 
        for j in range(0, n):
            x[j] = b[j]
            for k in range(0, j):
                x[j] = x[j] - a[j, k]*x[k]
            for k in range(j+1, n):
                x[j] = x[j] - a[j, k]*xprev[k]
            x[j] = x[j] / a[j, j]
 
        # Start på konvergenstest.
        i = 0
        for j in range(0, n):
            if abs(x[j] - xprev[j]) > abs(x[i] - xprev[i]):
                i = j
        if abs(x[i] - xprev[i]) < eps*abs(x[i]):
            done = True
        # Slutt på konvergenstest.
 
    return x

Si at vi kjører koden med

a = np.array([[1.0,   -0.25, -0.25, 0    ],
              [-0.25, 1,     0,     -0.25],
              [-0.25, 0,     1,     -0.25],
              [0,     -0.25, -0.25, 1    ]])
b = np.array([[50.0],
              [50.0],
              [25.0],
              [25.0]])
x0 = np.array([[1.0],
               [1.0],
               [1.0],
               [1.0]])
eps = 0.00001
 
foo(a, b, x0, eps)

Hvis vi antar uendelig maskinpresisjon (altså at flyttallsrepresentasjon ikke skaper trøbbel), kan du da være sikker på at koden terminerer (altså at done settes til True på et eller annet tidspunkt)? Hvorfor? Hva nærmer returverdien seg når eps gjøres liten? (Du trenger ikke å finne returverdien, bare forklare hvilken spesiell vektor den er.) Begrunn svarene! (Hint: Begrunnelsen avhenger av den angitte matrisen a).