Review

Last week we looked at a way to solve for unknown head with the following dirichlet boundary conditions. (If you haven’t seen the first post look here before continuing)

boundary_conditions <- c(8.04, 8.18, 8.36, 8.53,
                        7.68, 0,    0,    8.41,
                        7.19, 0,    0,    8.33,
                        6.82, 7.56, 7.99, 8.29)

A <- matrix(data = boundary_conditions, nrow = 4, ncol = 4, byrow = TRUE)

print(A, digits = 3)
##      [,1] [,2] [,3] [,4]
## [1,] 8.04 8.18 8.36 8.53
## [2,] 7.68 0.00 0.00 8.41
## [3,] 7.19 0.00 0.00 8.33
## [4,] 6.82 7.56 7.99 8.29

Variations on a theme

We were able to solve this problem with a method called Gauss-Seidel iteration. This time, we going to try a small variation on that method called successive over relaxation.

Again, we’ll create indices for our matrix and take some initial guesses at the unknown values (i.e. zeros from above)

B <- matrix(seq_along(A), nrow = 4, ncol = 4, byrow = TRUE)
B
##      [,1] [,2] [,3] [,4]
## [1,]    1    2    3    4
## [2,]    5    6    7    8
## [3,]    9   10   11   12
## [4,]   13   14   15   16
# get rid of the matrix boundaries and focus on the unknown values
B <- B[c(-1, -4), c(-1, -4)]
B
##      [,1] [,2]
## [1,]    6    7
## [2,]   10   11
A[c(B)] <- c(8, 8.5, 7, 8)

# Check that are initial conditions and guesses are good to go
print(A, digits = 3)
##      [,1] [,2] [,3] [,4]
## [1,] 8.04 8.18 8.36 8.53
## [2,] 7.68 8.00 8.50 8.41
## [3,] 7.19 7.00 8.00 8.33
## [4,] 6.82 7.56 7.99 8.29

Looks good. We’re ready to go.

Solving the problem

Again, we will solve this problem with iteration, but now we’re going to add an extra parameter called “omega” that will hopefully help us reach our goal faster. So let’s review our additional inputs to this algorithm:

This time, we’ll keep track of all of our residuals so that we can visualize our path to convergence. Again, we need to create an appropriately sized matrix to store these results.

# Only let the loop go 100 times
max_iteration <- 100
# Stop the loop if the change between old and new values is less than 0.001
convergence_criteria <- 1e-3

# a small matrix to store the differences between old and new values.
# We can keep them all in a single column. We'll start this just above 
# our convergence criteria, but they will get replaced inside the inner for loop

residual <- matrix(0.0011, length(B)*max_iteration, 1) 
# This matrix is as long as the max number of iterations
# times the number of unknown head values

# the maximum difference between old and new values will serve as our convergence criteria
err <- max(residual)

# the omega value (The heuristic is to set this somewhere between 1-2)
# let's start with 1.2
omega <- 1.2

# Counter to use for storing our residuals in the matrix we just created
a <- 1

Now we are ready to solve this problem

We will solve this with a loop within a loop.

for(t in 1:max_iteration) {                                           # Outer loop iterates over the matrix
                                                                      # t counts the outer loop iteration number
  if(err >= convergence_criteria) {                                   # Continue loop if we haven't converged
  
                                                            
  
  for(i in B) {                                                       # Inner loop that replaces guesses
    
    old_value <- A[i]                                                 # What was originally in the matrix
    if(t > 1){
      # Here is the SOR equation. It's just our 5-point solution plus an additional term 
      # that is proportional to our residual time our omega parameter.
      # Thus, extrapolation is greater if the previous residual was large
    A[i] <- ((A[i-1]+A[i-length(A[,2])]+A[i+1]+A[i+length(A[,2])])/4)+(omega*residual[a-4])
                                                                      # extrapolating 
    } else {
      # The 5-point solution. Just plain old Gauss-Seidel
    A[i] <- ((A[i-1]+A[i-length(A[,2])]+A[i+1]+A[i+length(A[,2])])/4)
    
    }  
    residual[a] <- A[i]-old_value                   # Difference between old and new
    a <- a + 1                                      # Increment our counter
  }
  
  err <- max(abs(residual[(t*4-3):(t*4)]))          # Set our criteria to the max residual
  
  } else if(err < convergence_criteria) {           # If we get below the criteria, 
    break                                           # break the loop early.
  }
}

if(err < convergence_criteria){
print(paste(t-1, "iterations to reach convergence criteria"))
} else {print("convergence criteria not met :( Try adjusting omega")}
## [1] "convergence criteria not met :( Try adjusting omega"

Results

Oh no! A frown face… I thought successive overrelaxation was supposed to help us. Let’s look at our final results.

print(A, digits = 3)
##      [,1]      [,2]      [,3] [,4]
## [1,] 8.04  8.18e+00  8.36e+00 8.53
## [2,] 7.68  1.24e+08 -1.08e+08 8.41
## [3,] 7.19 -1.08e+08 -3.83e+08 8.33
## [4,] 6.82  7.56e+00  7.99e+00 8.29

Clearly something went very wrong here. Let’s plot up the residuals to see what’s happening. I’m going to hide my plotting code for now because that isn’t the point of this post. Let’s just focus on what omega is doing here.

SOR_plot_omega_1.2

Interesting, it’s hard to say for sure because the scale is so huge, but it looks like we may have been close to convergence but then actually diverged.

Let’s rerun this with a different omega value and see what happens. Again, I’m hiding the code to reset our initial conditions and rerun the analysis so we can focus on results from changes in omega.

Let’s break from the heuristic and try a value below 1 to see what happens.

omega <- 0.8
## [1] "convergence criteria not met :( Try adjusting omega"

Frowny face again?! Let’s take a look.

SOR_plot_omega_0.8

OK. Now we can clearly see how SOR almost converged early on, but overshot the convergence criteria and began oscillating in a growing sinusoidal pattern.

Let’s go nuts and try a value of 0.25 to see what happens.

omega <- 0.25
## [1] "10 iterations to reach convergence criteria :)"

Hooray! We finally converged and managed to get a smiley face! However, SOR still took longer to converge than just plain old Gauss-Siedel iteration. Let’s look at this plot to diagnose what’s happening.

SOR_plot_omega_0.25

If we set omega to 0, the additional SOR term on the right hand side of our equation becomes 0 and we’re back Gauss-Seidel. Let’s do that and then put the plots on top of each other for comparison.

omega <- 0
## [1] "5 iterations to reach convergence criteria :)"
library(patchwork)

SOR_plot_omega_0.25/SOR_plot_omega_0.0

Questions

This clearly shows how SOR has a tendency to overshoot the convergence criteria. The consequence of this is SOR can create oscillations that take extra time to settle.

Some questions to think about.

  1. If SOR performs worse the Gauss-Seidel, why use it?
  2. What would have happened if our convergence criteria was less strict?
  3. What would have happened if we had worse initial guesses?
  4. In what kind of situation would SOR be beneficial instead of harmful?