R Notebook

Sistemas de Ecuaciones Lineales

LU<-function(A) {
n<-nrow(A)
l<-rep(0,n^2)
dim(l)<-c(n,n)
u<-l
for(k in 1:n) {
for(j in k:n) {
u[k,j]<-A[k,j]
if(k > 1) {
for (p in 1:(k-1)) {
u[k,j]<-u[k,j]-l[k,p]*u[p,j]
}
}
}
for (i in k:n) {
l[i,k]<-A[i,k]
if (k > 1) {
for (p in 1:(k-1)) {
l[i,k]<-l[i,k]-l[i,p]*u[p,k]
}
}
l[i,k]<-l[i,k]/u[k,k]
}
}
return(list(L=l,U=u))
}

Revisar esto para la matriz A del taller

A = matrix(c(4,1,2,1,3,2,1,5,7),nrow=3)
A

     [,1] [,2] [,3]
[1,]    4    1    1
[2,]    1    3    5
[3,]    2    2    7

LU(A)

$L
     [,1]      [,2] [,3]
[1,] 1.00 0.0000000    0
[2,] 0.25 1.0000000    0
[3,] 0.50 0.5454545    1

$U
     [,1] [,2]     [,3]
[1,]    4 1.00 1.000000
[2,]    0 2.75 4.750000
[3,]    0 0.00 3.909091

Metodo Directo de valores Característicos

Raices y vectores caracteristicos

Método directo para hallar el vector y valor propio mas grande.

a <- cbind(c(5, -4, 1, 0),c(-4,6,-4,1),c(1,-4,6,-4),c(0,1,-4,5))
# valor inicial
x<-c(0,1,0,0)
lamda1 <-0
lamda <-max(a)
#Inicio del proceso iterativo
for (i in 1:5) {
lamda1 <- lamda
y <- a%*%x
lamda <- max(y)
x <- y/lamda
print(lamda)
z <- cbind(x,y)
print(z)
}

[1] 6
           [,1] [,2]
[1,] -0.6666667   -4
[2,]  1.0000000    6
[3,] -0.6666667   -4
[4,]  0.1666667    1
[1] 11.5
           [,1]      [,2]
[1,] -0.6956522 -8.000000
[2,]  1.0000000 11.500000
[3,] -0.8115942 -9.333333
[4,]  0.3913043  4.500000
[1] 12.42029
           [,1]       [,2]
[1,] -0.6674446  -8.289855
[2,]  1.0000000  12.420290
[3,] -0.8961494 -11.130435
[4,]  0.4994166   6.202899
[1] 12.75379
           [,1]       [,2]
[1,] -0.6455627  -8.233372
[2,]  1.0000000  12.753792
[3,] -0.9441903 -12.042007
[4,]  0.5552608   7.081680
[1] 12.91427
           [,1]       [,2]
[1,] -0.6327885  -8.172004
[2,]  1.0000000  12.914273
[3,] -0.9703797 -12.531747
[4,]  0.5848618   7.553065

metodo Iterativo

[1] -2.0  0.1
[1] 0.5257311 0.1000000
[1] -0.8464782  0.1000000
[1] -0.0712627  0.0100000
[1] -0.04474848  0.01000000
[1] -0.0007983883  0.0001000000
[1] 6.082346e-06 1.000000e-06
[1] 2.585073e-09 1.000000e-09
[1] -3.935092e-15  1.000000e-15
$raices
[1] 6.2899451 2.2942804 0.4157746

$vectores
           [,1]       [,2]      [,3]
[1,]  0.8359921  0.5048961 0.2149353
[2,] -0.5391919  0.6830536 0.4926559
[3,]  0.1019277 -0.5277478 0.8432633

$estado
              [,1]          [,2]          [,3]
[1,]  6.289945e+00 -1.731734e-24 -3.627635e-39
[2,] -1.731734e-24  2.294280e+00  0.000000e+00
[3,] -3.627635e-39  0.000000e+00  4.157746e-01

Algunos funciones en R asociados a criterios de convergencia

a <- cbind(c(5, -4, 1, 0),c(-4,6,-4,1),c(1,-4,6,-4),c(0,1,-4,5))
espectral <- eigen(a)
vectores <- as.matrix(espectral$vectors)
valores <- as.matrix(espectral$values)
print(a)

     [,1] [,2] [,3] [,4]
[1,]    5   -4    1    0
[2,]   -4    6   -4    1
[3,]    1   -4    6   -4
[4,]    0    1   -4    5

print(valores)

          [,1]
[1,] 13.090170
[2,]  6.854102
[3,]  1.909830
[4,]  0.145898

print(vectores)

          [,1]      [,2]      [,3]     [,4]
[1,]  0.371748  0.601501 -0.601501 0.371748
[2,] -0.601501 -0.371748 -0.371748 0.601501
[3,]  0.601501 -0.371748  0.371748 0.601501
[4,] -0.371748  0.601501  0.601501 0.371748

Soluciones del Taller

2. Dada la siguiente matriz, utilice las funciones del paquete para descomponer la matriz \(A=L+D+U\) (Jacobi)

A = matrix(c(-8.1, -7, 6.123, -2, -1, 4,
-3, -1, 0, -1, -5, 0.6,
-1, 0.33, 6, 1/2), nrow=4, byrow=TRUE)
A

     [,1]  [,2]   [,3] [,4]
[1,] -8.1 -7.00  6.123 -2.0
[2,] -1.0  4.00 -3.000 -1.0
[3,]  0.0 -1.00 -5.000  0.6
[4,] -1.0  0.33  6.000  0.5

PUNTO 2A

require(Matrix)
require(pracma)
A = matrix(c(-8.1, -7, 6.123, -2, -1, 4,
             -3, -1, 0, -1, -5, 0.6,
             -1, 0.33, 6, 1/2), nrow=4, byrow=TRUE)
A

     [,1]  [,2]   [,3] [,4]
[1,] -8.1 -7.00  6.123 -2.0
[2,] -1.0  4.00 -3.000 -1.0
[3,]  0.0 -1.00 -5.000  0.6
[4,] -1.0  0.33  6.000  0.5

ludec = lu(A)
L <- ludec$L
L

          [,1]       [,2]      [,3] [,4]
[1,] 1.0000000  0.0000000  0.000000    0
[2,] 0.1234568  1.0000000  0.000000    0
[3,] 0.0000000 -0.2055838  1.000000    0
[4,] 0.1234568  0.2455076 -1.068263    1

U <-ludec$U
U

     [,1]      [,2]      [,3]       [,4]
[1,] -8.1 -7.000000  6.123000 -2.0000000
[2,]  0.0  4.864198 -3.755926 -0.7530864
[3,]  0.0  0.000000 -5.772157  0.4451777
[4,]  0.0  0.000000  0.000000  1.4073689

D <- diag(diag(A))
D

     [,1] [,2] [,3] [,4]
[1,] -8.1    0    0  0.0
[2,]  0.0    4    0  0.0
[3,]  0.0    0   -5  0.0
[4,]  0.0    0    0  0.5

A = L %*% U
A

     [,1]  [,2]   [,3] [,4]
[1,] -8.1 -7.00  6.123 -2.0
[2,] -1.0  4.00 -3.000 -1.0
[3,]  0.0 -1.00 -5.000  0.6
[4,] -1.0  0.33  6.000  0.5

2. Utilice la función itersolve(A, b, tol , method = “Gauss-Seidel”) y solucionar el sistema asociado a la matriz \(A\) con \(b=[1.45,3,5.12,-4]^{t}\) con una tolerancia de \(1e^-9\)

PUNTO 2B

A = matrix(c(-8.1, -7, 6.123, -2, -1, 4,
             -3, -1, 0, -1, -5, 0.6,
             -1, 0.33, 6, 1/2), nrow=4, byrow=TRUE)

b <- matrix(c(1.45,3,5.12,4.0), nrow = 4, ncol = 1, byrow = TRUE)

cat("Mediante m?todo de Gauss-Seidel\n")

itersolve(A, b, tol = 1e-9, method = "Gauss-Seidel")

Salida: Iteraci?n 992 -> Error relativo = 8.45536e+195 Iteraci?n 993 -> Error relativo = 1.330449e+196 Iteraci?n 994 -> Error relativo = 2.093458e+196 Iteraci?n 995 -> Error relativo = 3.294051e+196 Iteraci?n 996 -> Error relativo = 5.18318e+196 Iteraci?n 997 -> Error relativo = 8.155721e+196 Iteraci?n 998 -> Error relativo = 1.283301e+197 Iteraci?n 999 -> Error relativo = 2.01927e+197 Iteraci?n 1000 -> Error relativo = 3.177317e+197

N?mero de iteraciones realizado fue 1000

Soluciones:[1] -9.448624e+196 7.294493e+196 4.118582e+196 -7.313460e+197

3. Genere 5 iteraciones del método de Jacobi, calcular error relativo para cada iteracion

PUNTO 2C MODIFICANDO FUNCI?N ITERSOLVE

itersolve <- function(A, b, x0 = NULL, 
                      nmax = 1000, tol = .Machine$double.eps^(0.5),
                      method = c("Gauss-Seidel", "Jacobi", "Richardson")) {
  stopifnot(is.numeric(A), is.numeric(b))
  
  n <- nrow(A)
  if (ncol(A) != n)
    stop("Argument 'A' must be a square, positive definite matrix.")
  b <- c(b)
  if (length(b) != n)
    stop("Argument 'b' must have the length 'n = ncol(A) = nrow(A).")
  if (is.null(x0)) {
    x0 <- rep(0, n)
  } else {
    stopifnot(is.numeric(x0))
    x0 <- c(x0)
    if (length(x0) != n)
      stop("Argument 'x0' must have the length 'n=ncol(A)=nrow(A).")
  }
  
  method <- match.arg(method)
  
  if (method == "Jacobi") {
    L <- diag(diag(A))
    U <- eye(n)
    beta <- 1; alpha <- 1
  } else if (method == "Gauss-Seidel") {
    L <- tril(A)
    U <- eye(n)
    beta <- 1; alpha <- 1
  } else {  # method = "Richardson"
    L <- eye(n)
    U <- L
    beta <- 0
  }
  
  b <- as.matrix(b)
  x <- x0 <- as.matrix(x0)
  r <- b - A %*% x0
  r0 <- err <- norm(r, "f")
  
  iter <- 0
  while (err > tol && iter < nmax) {
    iter <- iter + 1
    z <- qr.solve(L, r)
    z <- qr.solve(U, z)
    if (beta == 0) alpha <- drop(t(z) %*% r/(t(z) %*% A %*% z))
    x <- x + alpha * z
    r <- b - A %*% x
    err <- norm(r, "f") / r0
    cat("Iteraci?n ",iter," -> Error relativo = ", err,"\n")
  }
  
  cat ("\nN?mero de iteraciones realizado fue ",iter)
  
  cat ("\n\nSoluciones:")
  print(c(x))
}
A = matrix(c(-8.1, -7, 6.123, -2, -1, 4,
             -3, -1, 0, -1, -5, 0.6,
             -1, 0.33, 6, 1/2), nrow=4, byrow=TRUE)
b <- matrix(c(1.45,3,5.12,4.0), nrow = 4, ncol = 1, byrow = TRUE)
cat("Mediante m?todo de Jacobi\n")

Mediante m?todo de Jacobi

itersolve(A, b, nmax = 5, tol = 1e-9, method = "Jacobi")

Iteraci?n  1  -> Error relativo =  3.943147 
Iteraci?n  2  -> Error relativo =  4.117929 
Iteraci?n  3  -> Error relativo =  4.585134 
Iteraci?n  4  -> Error relativo =  9.205016 
Iteraci?n  5  -> Error relativo =  9.634505 

N?mero de iteraciones realizado fue  5

Soluciones:[1]  3.177583 -6.415903 -3.440095 20.070249

PUNTO 2C IMPLEMENTANDO FUNCION DE NOVO

library( pracma)
A <- matrix(c(-8.1, -7, 6.123, -2, -1, 4,
              -3, -1, 0, -1, -5, 0.6,
              1, 0.33, 6, 1/2), nrow=4, byrow=TRUE)
print(A)

     [,1]  [,2]   [,3] [,4]
[1,] -8.1 -7.00  6.123 -2.0
[2,] -1.0  4.00 -3.000 -1.0
[3,]  0.0 -1.00 -5.000  0.6
[4,]  1.0  0.33  6.000  0.5

b<-c(8,15,1,-4) 
f1<-function(A,x,B,n){
  sup<-matrix(0,n,n)
  inf<-matrix(0,n,n)
  diag<-matrix(0,n,n)
  for (i in 0:n){
    for (j in 0:n){
      if(j>i){
        sup[i,j]=A[i,j]
      }
      if(j<i){
        inf[i,j]=A[i,j]
      }
      if(j==i){
        diag[i,j]=A[i,j]
      }
    }
  }
  print(" descomposicion LDU ")
  print("diagonal superior")
  print("  ")
  print(sup)
  print("  ")
  print("diagonal inferior")
  print("  ")
  print(inf)
  print("  ")
  print("diagonal ")
  print("  ")
  print(diag)
  f2(sup,inf,diag,x,B)
  
}
f2<-function(sup,inf,diag,x,B){
  n=0
  print("iteracion ")
  print(n)
  inv=solve(diag)
  err=1
  while( n<5){
    print("iteraci?n ")
    print(n)
    x1<-x
    x<-(inv*(B-(inf+sup)*x1))
    print(x)
    err=(x[1]-x1[1])/x[1]
    print("error: ")
    print(err)
    n=n+1
  }
  
  print("prueba final: ")
  print((inf+sup+diag)*x)
}
x<-c(1,2,3,1)
f1(A,x,b,4)

[1] " descomposicion LDU "
[1] "diagonal superior"
[1] "  "
     [,1] [,2]   [,3] [,4]
[1,]    0   -7  6.123 -2.0
[2,]    0    0 -3.000 -1.0
[3,]    0    0  0.000  0.6
[4,]    0    0  0.000  0.0
[1] "  "
[1] "diagonal inferior"
[1] "  "
     [,1]  [,2] [,3] [,4]
[1,]    0  0.00    0    0
[2,]   -1  0.00    0    0
[3,]    0 -1.00    0    0
[4,]    1  0.33    6    0
[1] "  "
[1] "diagonal "
[1] "  "
     [,1] [,2] [,3] [,4]
[1,] -8.1    0    0  0.0
[2,]  0.0    4    0  0.0
[3,]  0.0    0   -5  0.0
[4,]  0.0    0    0  0.5
[1] "iteracion "
[1] 0
[1] "iteraci?n "
[1] 0
           [,1] [,2] [,3] [,4]
[1,] -0.9876543 0.00  0.0    0
[2,]  0.0000000 3.75  0.0    0
[3,]  0.0000000 0.00 -0.2    0
[4,]  0.0000000 0.00  0.0   -8
[1] "error: "
[1] 2.0125
[1] "iteraci?n "
[1] 1
           [,1] [,2] [,3] [,4]
[1,] -0.9876543 0.00  0.0    0
[2,]  0.0000000 3.75  0.0    0
[3,]  0.0000000 0.00 -0.2    0
[4,]  0.0000000 0.00  0.0   -8
[1] "error: "
[1] 0
[1] "iteraci?n "
[1] 2
           [,1] [,2] [,3] [,4]
[1,] -0.9876543 0.00  0.0    0
[2,]  0.0000000 3.75  0.0    0
[3,]  0.0000000 0.00 -0.2    0
[4,]  0.0000000 0.00  0.0   -8
[1] "error: "
[1] 0
[1] "iteraci?n "
[1] 3
           [,1] [,2] [,3] [,4]
[1,] -0.9876543 0.00  0.0    0
[2,]  0.0000000 3.75  0.0    0
[3,]  0.0000000 0.00 -0.2    0
[4,]  0.0000000 0.00  0.0   -8
[1] "error: "
[1] 0
[1] "iteraci?n "
[1] 4
           [,1] [,2] [,3] [,4]
[1,] -0.9876543 0.00  0.0    0
[2,]  0.0000000 3.75  0.0    0
[3,]  0.0000000 0.00 -0.2    0
[4,]  0.0000000 0.00  0.0   -8
[1] "error: "
[1] 0
[1] "prueba final: "
     [,1] [,2] [,3] [,4]
[1,]    8    0    0    0
[2,]    0   15    0    0
[3,]    0    0    1    0
[4,]    0    0    0   -4

4. Cree una función que cuente el número de multiplicaciones en el método directo de Gauss Jordan, para resolver un sistema de \(n\) ecuaciones y pruebelo para \(n=5\)

PUNTO 4 EN R

n<-5
c<-matrix(0,n,n+1)
for (i in 0:n){
  for (j in 0:(n+1)){
    c[i,j]=floor(runif(1, min=0, max=101))
  }
  
}
contador<-function(matriz){
  tot<-0
  for (i in 1:nrow(matriz)){
    p<-matriz[i,i]
    for (j in 1:ncol(matriz)){
      matriz[i,j]<-matriz[i,j]/p
      tot<-tot+1
    }
    for (j in 1:nrow(matriz)){
      if(j!=i){
        p<-matriz[j,i]
        for (k in 1:ncol(matriz)){
          matriz[j,k]<-matriz[j,k]-p*matriz[i,k]
          tot<-tot+1
        }
        
      }
    }
  }
  print(matriz)
  print("N?mero de operaciones: ")
  print(tot)
}
print(c)

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]   16   42   77   64   56   75
[2,]   17   22    8   21   81   63
[3,]   31   60   22   75   94   89
[4,]    0   50   75   96   83   88
[5,]   33   41   29   40   26   24

contador(c)

     [,1] [,2] [,3] [,4] [,5]        [,6]
[1,]    1    0    0    0    0 -38.7561727
[2,]    0    1    0    0    0  63.3471474
[3,]    0    0    1    0    0   3.2265324
[4,]    0    0    0    1    0 -34.9956543
[5,]    0    0    0    0    1   0.4606688
[1] "N?mero de operaciones: "
[1] 150

5. Dado el siguiente sistema:

   \(2x-z=1\)
   \(\beta\)x+2y-z=2
   \(-x+y+\)\(z=1\)

1. Encuentre el valor de \(\alpha\) y \(\beta\) para asegura la convergencia por el m??todo de Jacobi

library(pracma)
beta <- 0
alpha<-0
A = matrix(c(2, 0, 1, beta,2 , -1,
-1, 1, alpha), nrow=3, byrow=TRUE)
A

     [,1] [,2] [,3]
[1,]    2    0    1
[2,]    0    2   -1
[3,]   -1    1    0

B = matrix (c(1,2,1),nrow=3, byrow=TRUE)
B

     [,1]
[1,]    1
[2,]    2
[3,]    1

beta<-0
alpha<-3