Let’s take a look at the Wang and Anderson (1982) example that uses a 5-point solution to solve for unknown groundwater head.
Here are the dirichlet boundary conditions to the problem. For more on boundary conditions related to ground water modeling, see the ModelMuse help page.
The zeros in the matrix below represent the unknown head, while the boundaries show constant head boundaries (i.e. dirichlet boundaries).
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.29Now let’s 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 11A[c(B)] <- c(8, 8.5, 7, 8)Now we have all of our initial conditions set up and ready 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.29We are going to solve this problem with iteration, but that means we need to decide a few things first.
The last thing we want is to accidentally create an infinite loop, so these conditions give us both control and safe guards.
We can also keep track of our progress towards convergence and see how long it takes to arrive at our final solution. We’ll refer to the difference between an old and new value as a residual. We don’t want to grow anything with in a loop (computationally inefficient) so we need to create an appropriately sized matrix to store these residual 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'll start this with values of 1,
# but they will get replaced inside the inner for loop
residual <- matrix(1, nrow = dim(B)[1], ncol = dim(B)[2]) # This matrix has the same dimensions as our matrix of unknowns
# the maximum difference between old and new values will serve as our convergence criteria
err <- max(residual)We will solve this with a loop within a loop.
for(t in 0: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
a <- 1 # Counter to use for storing our residuals
for(i in B) { # Inner loop that replaces guesses
old_value <- A[i] # What was originally in the matrix
A[i] <- (A[i-1]+A[i-length(A[,2])]+A[i+1]+A[i+length(A[,2])])/4 # The 5-point solution
residual[[a]] <- abs(A[i]-old_value) # Difference between old and new
a <- a + 1 # Increment our counter
}
err <- max(residual) # Set our criteria to the max residual
} else if(err < convergence_criteria) { # If we get below the criteria,
break # break the loop early.
}
}
print(paste(t, "iterations to reach convergence criteria"))## [1] "5 iterations to reach convergence criteria"Great, let’s look at our final results.
print(A, digits = 3)## [,1] [,2] [,3] [,4]
## [1,] 8.04 8.18 8.36 8.53
## [2,] 7.68 7.93 8.19 8.41
## [3,] 7.19 7.68 8.05 8.33
## [4,] 6.82 7.56 7.99 8.29