Programmering

Forlæns substitution - matrix i R

11. september kl. 17:04 af louisesørensen2

Hej derude,

Jeg skal skrive et program i R som løser et lineært ligningssystemt Lx=b ved forlæns substitution, hvor L er en nedre triangulær matrix. Input er matricen L og vektoren b, og output er løsningen x

vores bog snakker om at x_i=\frac{b_i-\sum_{j=1}^{i-1}a_{ij}x_j}{a_{ii}}, men jeg har simpelthen ikke nogen mulig idé hvordan jeg skal gribe den her opgave an... 

o.b.s. jeg er rigtig dårlig til at skrive programmer...


Brugbart svar (0)

Svar #1
11. september kl. 20:49 af norm

Skriv en funktion, som tager matricen L og vektoren b som input. Du kan starte med at definere en vektor fx x = b, og byt så det første element i vektoren ud med x[1]. Så kan du lave et for-loop, som beregner de næste x[i] og substituerer dem ind i den midlertidige x-vektor, du lavede til at starte med. Du skal altså implementere den viste formel i for-loopet.


Brugbart svar (0)

Svar #2
11. september kl. 20:51 af norm

Og jeg vil anbefale dig at lave være med at skrive et nested-loop (altså et loop i loopet), da det er svært at læse, fylder mere og kan gøre programmet langsommere, hvis matricerne bliver store. I stedet for kan du bruge det faktum, at du kan gange vektorer sammen elementvist i R. Du kan definere indekseringen i loopet til 2:nrow(L). Så kan du printe hver række af matricen ud ved at skrive L[i,] i loopet.


Svar #3
13. september kl. 16:56 af louisesørensen2

Hmm...

jeg har svært ved at se hvordan jeg kommer i gang, norm.

Jeg tænker

Lsolve<-function(L,b){

rækker<-nrow(L)

p<-1

for(i in 1:rækker){

}

}

her stopper den så for mig, fordi jeg vil jo gerne ha' at den laver guess-elimination på 1. række (a_11) også ved at der er 0'er i alle andre indgange. Hvordan fortæller jeg programmet det?

mvh


Brugbart svar (0)

Svar #4
13. september kl. 17:04 af norm

Nej, det er ikke Gausselimination, du skal lave. Det er forlæns substitution. Tilfældet er jo, at matricen er på nedre triangulær form, og derfor vil du undgå at lave Gausselimination, da Gausselimination forøger risikoen for ciffertab ved subtraktion og i øvrigt gør programmet langsommere. Du har jo et ligningssystem på formen:

a[1,1]*x1=b1

a[2,1]*x1+a[2,2]*x2=b2

Osv. Du starter ud med at finde x1=b1/a[1,1]. Det indsætter du i ligning nummer 2, og så isolerer du x2, og sådan fortsætter det. Der er selvfølgelig flere måder at skrive koden på, men et forslag er:

LTS <- function(L,b) {
  b[1] <- b[1]/L[1,1]
  
  for (i in 2:nrow(L)) {
     r <- L[i,]
     b[i] <- (b[i]-sum(r[-i]*b[-i]))/L[i,i]
    
  }
  return(b)
  
}


Svar #5
13. september kl. 17:10 af louisesørensen2

Alright, så du laver den rekursiv. Hvis det er ok, kan jeg så lige spørge hvad minus'et i indekset f.eks. r[-i] skal betyde?

Og i dit for loop siger du 2 fordi at du har "manuelt" lavet den første række?

hvis man vil have den til at return x vil du så bare sætte x = b[i]? og return(x)? 


Brugbart svar (0)

Svar #6
13. september kl. 17:23 af norm

Den er faktisk ikke rekursivt defineret, men iterativ. Men r[i] betyder, at du beder den om vektoren r minus det i'te element dvs. hvis R=[1,2,3], så er fx R[-2]=[1,3]. Og ja, for-loopet starter i 2, fordi man laver n-1 substitutioner, hvor n er antallet af rækker/søjler. Den returnerer faktisk x, men for at simplificere koden har jeg valgt bare at lade funktionen omskrive b frem for at definere ny vektor x, men det er selvfølgelig ikke så nemt at se.


Svar #7
13. september kl. 17:31 af louisesørensen2

Jeg kan simpelthen ikke få koden til at virke hmm.. den siger "Fejl i b[i] <- (b[i] - sum(r[-i] * b[-i]))/L[i, i] : 
  erstatning har længden nul". Har prøvet noget forskelligt.


Svar #8
13. september kl. 17:44 af louisesørensen2

Så jeg prøvede først bare at lave et program som kan løse nedre triangulær matricer som er dimension 3:

L3solve<-function(L,b){
  b[1]<-b[1]/L[1,1]
  b[2]<-(b[2]-L[2,1]*b[1])/L[2,2]
  b[3]<-(b[3]-L[3,1]*b[1]-L[3,2]*b[2])/L[3,3]
return(b)
  }

Til min store sucess virkede det (wuhu)

så prøvede jeg at kopierer programmet til en generaliserede program:

Lsolve<-function(L,b){
  b[1]<-b[1]/L[1,1]
  ræk<-nrow(L)
  for(i in 2:ræk){
    b[i]<-(b[i]-sum(L[i,i-1]*b[i-i]))/L[i,i]
  }
  return(b)
}

men det nægter simpelthen at virke... har du nogle forslag?


Brugbart svar (0)

Svar #9
13. september kl. 19:30 af norm

Det var mærkeligt. Mit program kører fint for mig. Er du sikker på, du har defineret matricen og vektoren korrekt? Når du kalder funktionen, skal du huske at give matricen som det første argument og vektoren som det næste argument.


Svar #10
13. september kl. 20:11 af louisesørensen2

Alt er defineret korrekt, har triple-tjekket.

Jeg kører det her program:

Lsolve<-function(L,b){
  b[1]<-b[1]/L[1,1]
  ræk<-nrow(L)
  for(i in 2:ræk){
    b[i]<-(b[i]-sum(L[i,]*b[i-1]))/L[i,i]
  }
  return(b)
}

og har følgende matrix og vektor:

A=\begin{bmatrix} 4 & 0 & 0\\ 2 & 2 &0 \\ 3& 3 & 2 \end{bmatrix}, \: \: \: c=\begin{pmatrix} 2\\ 1\\ 3 \end{pmatrix}

og når jeg kører programmet får jeg x=\begin{pmatrix} 0.5\\ 0.5\\ 1.5 \end{pmatrix}    hvilket er forskelligt fra når jeg kører R's eget forwardsolver funktion: forwardsolve(A,c)=\begin{pmatrix} 0.5\\ 0\\ 0.75 \end{pmatrix} 

jeg må gøre noget forkert? Kan simpelthen ikke forstå det.

Jeg lavede en funktion på baggrund af x_i=\frac{b_i-\sum_{j=1}^{i-1}a_{ij}x_j}{a_{ii}} som fungerer for de matricer der har dimension 3:

L3solve <- function(L,b) {
  b[1] <- b[1] / L[1,1]
  b[2] <- (b[2] - L[2,1]*b[1]) / L[2,2]
  b[3] <- (b[3] - L[3,1]*b[1] - L[3,2]*b[2]) / L[3,3]
  return(b)

og den virker, men jeg kan ikke forstå hvorfor generaliseringen ikke gider at virke optimalt.


Svar #11
13. september kl. 20:26 af louisesørensen2

Jeg kunne godt forestille mig at det er noget inden i summen som ikke gør hvad den skal. Betragter vi igen min hjemmestrikket Lsolve funktion for 3x3 matricer:

L3solve <- function(L,b) {
  b[1] <- b[1] / L[1,1]
  b[2] <- (b[2] - L[2,1]*b[1]) / L[2,2]
  b[3] <- (b[3] - L[3,1]*b[1] - L[3,2]*b[2]) / L[3,3]
  return(b)

ser vi jo netop at L[i,i-1]*b[1] skal blive til L[i,i-2]*b[2] osv. afhængig af hvilken række vi er på i matricen.

Jeg er kan ikke se om programmet gør det? 


Brugbart svar (0)

Svar #12
13. september kl. 21:07 af norm

Det her er fra din kode:

b[i]<-(b[i]-sum(L[i,]*b[i-1]))/L[i,i]

Problemet er, at du udtager hele rækken L[i,] og så ganger du kun ét element b[i-1] på rækken. Det fungerer ikke, fordi du skal ikke medtage det i'te element i summen. Desuden så skal du gange alle elementerne på nær det i'te i søjlevektoren b på rækkevektoren L[i,]. Læg mærke til, at hvis du fjerner det i'te element fra den i'te række, så får du en rækkevektor, hvor elementerne fra og med i+1 er 0'er. Grunden til at denne implementering er ok er, at det ikke giver anledning til afrundingsfejl i floating point systemet, fordi 0 jo er repræsenteret perfekt.


Skriv et svar til: Forlæns substitution - matrix i R

Du skal være logget ind, for at skrive et svar til dette spørgsmål. Klik her for at logge ind.
Har du ikke en bruger på Studieportalen.dk? Klik her for at oprette en bruger.