Forlæns substitution - matrix i R

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.

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.

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

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)
  
}

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)? 

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.

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.

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?

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.

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:

og når jeg kører programmet får jeg     hvilket er forskelligt fra når jeg kører R's eget forwardsolver funktion:  

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

Jeg lavede en funktion på baggrund af  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.

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? 

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.