R question solve

Linear Equation Solve

1) Write a function to solve a system of 2 linear equations in 2 variables (2x2 system).

cat("The equations are:\n a1x + b1y = c1\n a2x + b2y = c2\n\n")

The equations are:
 a1x + b1y = c1
 a2x + b2y = c2

a1 <- as.numeric(readline(prompt = "a1 = "))
4
b1 <- as.numeric(readline(prompt = "b1 = "))
3
c1 <- as.numeric(readline(prompt = "c1 = "))
1

a2 <- as.numeric(readline(prompt = "a2 = "))
5
b2 <- as.numeric(readline(prompt = "b2 = "))
7
c2 <- as.numeric(readline(prompt = "c2 = "))
8


values <- c(a1,b1,c1,a2,b2,c2)
eqn_2by2_solver <- function(a1, b1, c1, a2, b2, c2){
  # Solving the system
  denominator <- (a1*b2 - a2*b1)
  
  if (denominator == 0) {
    cat("No unique solution exists.\n")
  } else {
    y <- (a1*c2 - a2*c1) / denominator
    x <- (c1 - b1*y) / a1
    cat(sprintf("Solution:\nx = %.4f \ny = %.4f\n", x, y))
  }
}
eqn_2by2_solver(a1,b1,c1,a2,b2,c2)

Solution:
x = -1.3077 
y = 2.0769

do.call(eqn_2by2_solver, as.list(values))

Solution:
x = -1.3077 
y = 2.0769

2) Write a function to solve a 3x3 system using Cramer’s rule.

cramer3by3 <- function(eqn1,eqn2,eqn3){
  D <- rbind(
    eqn1[1:3],
    eqn2[1:3],
    eqn3[1:3]
  )
  b <- c(eqn1[4],eqn2[4],eqn3[4])
  vars <- c("x","y","z")
  for (i in 1:3){
    E <- D
    E[,i] <- b
    cat(vars[i],"=",det(E)/det(D),"\n")
  }
  
}

eqn1 <- c(a1=2, b1=3, c1=-1, d1=8)
eqn2 <- c(a2=1, b2=-1, c2=4, d2=11)
eqn3 <- c(a3=3, b3=1, c3=2, d3=7)

cramer3by3(eqn1,eqn2,eqn3)

x = -3.428571 
y = 6.714286 
z = 5.285714 

3) Write a function to find the determinant of a 2x2 matrix.

det2by2 <- function(A){
  (A[1,1]*A[2,2])-(A[1,2]*A[2,1])
}
A <- matrix(c(2,9,4,5), nrow=2)
det2by2(A)

[1] -26

4) Write a function to find the determinant of a 3x3 matrix

det3by3 <- function(A){
  result <- A[1,1]*(A[2,2]*A[3,3]-A[3,2]*A[2,3])-A[1,2]*(A[2,1]*A[3,3]-A[3,1]*A[2,3])+A[1,3]*(A[2,1]*A[3,2]-A[3,1]*A[2,2])
  return(result)
}

A=matrix(c(1:8,11), nrow=3)
det3by3(A)

[1] -6

A

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6   11

5) Write a function to find the inverse of a 2x2 matrix using the determinant and adjugate.

inv2by2 <- function(A){
  det_A <- (A[1,1]*A[2,2]-A[2,1]*A[1,2])
  if (det_A!=0){
  inv_A <- matrix(c(A[2,2], -A[1,2], -A[2,1], A[1,1]), nrow=2)/det_A
  return(inv_A)
  }
  else{
    return("Inverse not possible")
  }
}

A <- matrix(c(1, 2, 2, 3), ncol = 2, byrow = TRUE)
inv2by2(A)

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

6) Write a function to solve a 2x2 system by finding the inverse matrix and multiplying.

solver2by2 <- function(A, b){
  sol <- inv2by2(A) %*% b
  cat("x = ",sol[1,],"\ny = ", sol[2,])
}

A <- matrix(c(2,3,4,-1), byrow=TRUE, nrow=2)

b <- matrix(c(8,2), nrow=2)
b

     [,1]
[1,]    8
[2,]    2

solver2by2(A,b)

x =  1.142857 
y =  1.428571

7) Write a function to check if a 2x2 system has a unique solution, no solution, or infinitely many solutions.

For the system
a₁x + b₁y = c₁
a₂x + b₂y = c₂

Compute
D = a₁b₂ − a₂b₁

• If D ≠ 0 → unique solution
• If D = 0 and a₁/a₂ = b₁/b₂ = c₁/c₂ → infinitely many solutions
• If D = 0 and a₁/a₂ = b₁/b₂ ≠ c₁/c₂ → no solution

sol_type_2by2 <- function(A,b){
  D <- A[1,1]*A[2,2]-A[2,1]*A[1,2]
  if (D!=0){
    return ("Unique")
  }
  else{
    if (A[1,1]/A[2,1]==A[1,2]*A[2,2] & b[1,]/b[2,]==A[1,1]/A[2,1]){
      return("Infinitely many")
    } else{
      return("No Solution")
    }
  }
}

A <- matrix(c(2,3,4,-1), byrow=TRUE, nrow=2)
b <- matrix(c(8,2), nrow=2)
sol_type_2by2(A,b)

[1] "Unique"

8) Check the nature of roots of a quadratic equation

quadric_root_nature <- function(a,b,c){
  D <- b^2 - 4*a*c 
    if (D > 0) {
    cat("Two distinct real roots")
  } else if (D == 0) {
    cat("Two equal real roots")
  } else {
    cat("No real roots")
  }
}

a <- 1
b <- -3
c <- 2

quadric_root_nature(a,b,c)

Two distinct real roots

9) Write a function to check if a 3x3 system has a unique solution (by checking determinant).

sol_type_3by3 <- function(A,b){
  if (det(A)!=0){
    cat("Unique Solution")
  } else{
    cat("Infinitely many or No solution")
  }
}
sol_type_3by3(A,b)

Unique Solution

10) Write a function to find the trace of a square matrix (2x2, 3x3).

trace_finder <- function(A){
  if (nrow(A)==ncol(A)){
    t <- 0
    for (i in 1:nrow(A)){
      t <- t+A[i,i]
    }
    cat("Trace:",t)
  }
  else {
    cat("Give Square Matrix")
  }
}

A <- matrix(1:9, nrow = 3, ncol = 3)
trace_finder(A)

Trace: 15

matrix is singular when determinant = 0

---

Linear Algebra for Coders: Essential Concepts and Algorithms

1. Row Echelon Form (REF)

Definition

A matrix is in Row Echelon Form when:

• All nonzero rows are above rows of all zeros
• Each leading coefficient (pivot) of a row is in a column to the right of the leading coefficient of the row above it
• All entries in a column below a leading coefficient are zeros

Example

[ 1  2  3  4 ]
[ 0  1  5  6 ]   ✓ REF
[ 0  0  0  1 ]

[ 0  1  2 ]      ✗ Not REF (zero row above non-zero)
[ 1  0  3 ]

Algorithm to achieve REF

row_echelon_form <- function(mat) {
  rows <- nrow(mat)
  cols <- ncol(mat)
  pivot_row <- 1
  
  for (col in 1:cols) {
    # Find pivot in current column
    pivot <- find_pivot(mat, pivot_row, col)
    if (is.null(pivot)) {
      next
    }
    
    # Swap pivot row with current row
    mat <- swap_rows(mat, pivot, pivot_row)
    
    # Make all entries below pivot zero
    for (row in (pivot_row + 1):rows) {
      factor <- mat[row, col] / mat[pivot_row, col]
      mat[row, ] <- mat[row, ] - factor * mat[pivot_row, ]
    }
    
    pivot_row <- pivot_row + 1
    if (pivot_row > rows) {
      break
    }
  }
  return(mat)
}

---

2. Forward Elimination

Purpose

Convert a system of linear equations to upper triangular form.

Process

1. Start with first equation
2. Use it to eliminate first variable from subsequent equations
3. Move to next equation, repeat

Example System

2x +  y +  z =  8
4x + 3y + 2z = 18
-2x +  y +  z =  2

After Forward Elimination

2x +  y +  z = 8
0x +  y + 0z = 2  (Second equation minus 2×first)
0x + 2y + 2z = 10 (Third equation plus first)

---

3. Gaussian Elimination

Complete Algorithm

1. Forward Elimination: Convert to upper triangular matrix
2. Back Substitution: Solve from bottom to top

Implementation Steps

gaussian_elimination <- function(A, b) {
  n <- nrow(A)
  
  # Augment matrix [A|b]
  Ab <- cbind(A, b)
  
  # Forward elimination
  for (i in 1:n) {
    # Partial pivoting for stability
    max_row <- which.max(abs(Ab[i:n, i])) + i - 1
    if (i != max_row) {
      Ab[c(i, max_row), ] <- Ab[c(max_row, i), ]
    }
    
    # Eliminate below pivot
    for (j in (i + 1):n) {
      factor <- Ab[j, i] / Ab[i, i]
      Ab[j, ] <- Ab[j, ] - factor * Ab[i, ]
    }
  }
  
  # Back substitution
  x <- numeric(n)
  for (i in n:1) {
    x[i] <- Ab[i, n + 1] / Ab[i, i]
    for (j in (i - 1):1) {
      Ab[j, n + 1] <- Ab[j, n + 1] - Ab[j, i] * x[i]
    }
  }
  
  return(x)
}

---

---

4. Rank Check

Definition

The rank of a matrix is the maximum number of linearly independent rows (or columns).

Properties

• rank(A) ≤ min(m, n) for m×n matrix
• Full rank: rank = min(m, n)
• Singular matrix: rank < n for n×n matrix

Algorithm to Compute Rank

matrix_rank <- function(mat, augmented = FALSE) {
  ref_mat <- row_echelon_form(mat)
  rank <- 0
  cols <- ifelse(augmented, ncol(ref_mat) - 1, ncol(ref_mat))
  
  for (i in 1:nrow(ref_mat)) {
    if (any(abs(ref_mat[i, 1:cols]) > 1e-10)) {
      rank <- rank + 1
    }
  }
  return(rank)
}

Significance

• Unique solution: rank(A) = rank([A|b]) = n
• Infinite solutions: rank(A) = rank([A|b]) < n
• No solution: rank(A) < rank([A|b])

---

5. Augmented Matrix

Definition

Combines coefficient matrix A and constant vector b into [A|b]

Purpose

Allows simultaneous operations on both sides during elimination.

Example

System:

2x + 3y = 8
4x - y = 2

Augmented Matrix:

[ 2  3 | 8 ]
[ 4 -1 | 2 ]

Implementation

augment_matrix <- function(A, b) {
  cbind(A, b)
}

---

6. Solution Types for Linear Systems

6.1 Unique Solution

• Condition: rank(A) = rank([A|b]) = n
• Determinant: det(A) ≠ 0
• Interpretation: System has exactly one solution
• Example: Intersection of two non-parallel lines

6.2 Inconsistent System (No Solution)

• Condition: rank(A) < rank([A|b])
• Interpretation: Equations contradict each other
• Example: Parallel lines that don’t intersect

x + y = 3
x + y = 5  # Contradiction!

6.3 Infinite Solutions

• Condition: rank(A) = rank([A|b]) < n
• Free variables: n - rank(A)
• Interpretation: System has infinitely many solutions
• Example: Two equations representing the same line

---

7. Gauss-Jordan Elimination

Difference from Gaussian Elimination

• Creates Reduced Row Echelon Form (RREF)
• Both forward elimination AND backward elimination
• Pivot elements become 1
• All elements above and below pivots become 0

RREF Conditions

1. In REF
2. All pivot elements = 1
3. Each pivot is the only non-zero entry in its column

Algorithm

gauss_jordan <- function(A, b) {
  n <- nrow(A)
  Ab <- cbind(A, b)
  
  for (i in 1:n) {
    pivot <- Ab[i, i]
    for (j in i:(n + 1)) {
      Ab[i, j] <- Ab[i, j] / pivot
    }
    
    for (k in 1:n) {
      if (k != i) {
        factor <- Ab[k, i]
        for (j in i:(n + 1)) {
          Ab[k, j] <- Ab[k, j] - factor * Ab[i, j]
        }
      }
    }
  }
  
  return(Ab[, n + 1])
}

---

---

🎯 Gaussian Solver: Complete Visual Workflow

Gaussian Workflow

📋 Problem Overview

Solve Ax = b using Gaussian elimination, determining solution type without built-in functions.

# System to solve:
# 2x₁ + x₂ = 5
# x₁ + 3x₂ = 6

A <- matrix(c(2, 1, 1, 3), nrow = 2, byrow = TRUE)
b <- c(5, 6)

---

🏗️ System Representation

Original System:
┌               ┐
│  2   1  |  5  │  Equation 1: 2x₁ + x₂ = 5
│  1   3  |  6  │  Equation 2: x₁ + 3x₂ = 6
└               ┘

---

---

🧮 Final Implementation with Comments

gaussian_solver <- function(A, b) {
  n <- nrow(A)
  
  # Create augmented matrix manually (without cbind)
  aug <- matrix(0, nrow = n, ncol = n + 1)
  for (i in 1:n) {
    for (j in 1:n) {
      aug[i, j] <- A[i, j]
    }
    aug[i, n + 1] <- b[i]
  }
  
  # Forward elimination (no pivoting for simplicity)
  for (i in 1:n) {
    if (abs(aug[i, i]) < 1e-10) {
      next
    }
    
    for (j in (i + 1):n) {
      factor <- aug[j, i] / aug[i, i]
      for (k in i:(n + 1)) {
        aug[j, k] <- aug[j, k] - factor * aug[i, k]
      }
    }
  }
  
  # Check solution type
  has_zero_row <- FALSE
  has_inconsistent <- FALSE
  
  for (i in 1:n) {
    all_zero <- all(abs(aug[i, 1:n]) < 1e-10)
    if (all_zero) {
      has_zero_row <- TRUE
      if (abs(aug[i, n + 1]) > 1e-10) {
        has_inconsistent <- TRUE
      }
    }
  }
  
  if (has_inconsistent) {
    return(list(status = "No solution", solution = NULL))
  }
  
  if (has_zero_row) {
    return(list(status = "Infinite solutions", solution = NULL))
  }
  
  # Back substitution
  x <- numeric(n)
  for (i in n:1) {
    sum_val <- aug[i, n + 1]
    for (j in (i + 1):n) {
      sum_val <- sum_val - aug[i, j] * x[j]
    }
    x[i] <- sum_val / aug[i, i]
  }
  
  names(x) <- paste0("x", 1:n)
  return(list(status = "Unique solution", solution = x))
}

# Test the function
A <- matrix(c(2, 1, 1, 3), nrow = 2, byrow = TRUE)
b <- c(5, 6)
result <- gaussian_solver(A, b)
result$status
result$solution

---

🎯 Key Learning Points

1. Augmented Matrix: Combine A and b for simultaneous operations

2. Partial Pivoting: Prevents division by zero and improves numerical stability

3. Forward Elimination: Creates upper triangular form

4. Solution Analysis:

   • Inconsistent if [0 0 ... 0 | c] with c ≠ 0
   • Infinite if rank(A) < n
   • Unique if rank(A) = n and consistent

5. Back Substitution: Solve from bottom to top

---