Coding Assignment 4: System of Equations

Problem 1

Analytical Solution

First, the values given in the problem were rewritten in matrix notation and as a system of equations modeling the production run. The analytical solution was found using the backslash solver which is a built-in MATLAB command for solving a system of equations. Given the weekly total amount of copper, zinc, and glass available, 120 transistors, 100 resistors, and 90 computer chips will be manufactured this week. The problem setup and solution is shown below.

Numerical Solution

In MATLAB, Gauss Elimination with partial pivoting was employed to solve the system of equations. The results of the numerical method are shown below.

Sample Calculation Using Gauss Elimination With Partial Pivoting

Forward elimination for the first column and the subsequent back substitution is shown. In a complete calculation, forward elimination would be completed for each row before using back substitution so that the solutions can be easily found.

Sample Calculation Using The Gauss-Seidel Method (Extra Credit)

In addition to using Gauss Elimination with partial pivoting, the system was solved by hand using the Gauss-Seidel Method. This iterative approach starts with an initial guess and then systematically improves the estimate of the solution. After each iteration, the error between the current estimate and the last estimate is calculated. In real-world scenarios, there would be an error criteria that the solution must satisfy but in this problem, four iterations of the Gauss-Seidel Method are calculated. After just four iterations, the relative percent error of the number of transistors, resistors, and computer chips is less than 10%. Carrying out more iterations will increase the accuracy.

Error Analysis

True error between the analytical solution and the numerical solution was calculated in MATLAB using the following formula.

To confirm the true error between the analytical and numerical solution, the "format long" command was used. After using more significant digits, the true error remained the same.

Next, the system was checked for ill-conditioning. Using the inverse of [A], three methods for determining whether the system is ill-conditioned were performed. The first method requires that matrix [A] be scaled such that the highest value in each row is 1. If the inverse of this scaled matrix has values that are much greater than 1, then the system is ill-conditioned. Secondly, if the product of the inverse of the original [A] matrix and A is not the identity matrix, the system is ill-conditioned. The third method for checking ill-conditioning is taking the inverse of the inverse of [A]. This should yield the original [A]. If not, the system is ill-conditioned. All three methods do not indicate that the system is ill-conditioned.

Lastly, the condition number of the system was calculated using MATLAB. A condition number is a measure of the change in the output of a system given proportional to the input. Greater condition numbers indicate that more precision in the system is lost and the system is likely ill-conditioned. In this system, the small condition number of 8.5263 is consistent with the previous conditioning methods, indicating that the system is well-conditioned.

Discussion

Gauss elimination with partial pivoting accurately solved the system of equations to provide the correct number of transistors, resistors, and computer chips. The absolute true error between the solution found using Gauss elimination with partial pivoting and MATLAB's backslash solver was zero. The high accuracy can be attributed to partial pivoting. Partial pivoting avoids division by zero, minimizes round-off error, and improves ill-conditioning if present. Gauss elimination without partial pivoting should not be employed to solve a complex system of equations due to these advantages.

The condition number of a system should always be calculated to provide an indication of the conditioning. An ill-conditioned system of linear equations occurs when one or more of the equations have similar or the same slope. In this case, the intersection of the lines, or the solution, is difficult or impossible to pinpoint. On the contrary, well-conditioned systems have equations with dissimilar slopes. After using a numerical approach for solving a system of equations, finding the condition number is an important tool that can provide an explanation to the inaccuracy of the solution.

Pseudocode

Set Parameters A, b, nAnalytical Solutionbackslash solverr=0
DOFOR col = 1, n-1
DOIF A_col,col == 0row = col
DOFOR row = row +1, nENDDOcol1 = A_col,lengthb1 = b_colA_col,length = A_row,lengthb_col = b_rowA_row,length = col1b_row = b1ENDDO
DOFOR j = 1+r,n-1C = A_j+1,col / A_col,colA_j+1,length = A_j+1,length - C * A_col,lengthb_j+1 = b_j+1 - C * b_colENDDOr = r +1ENDDO
Substitutex_n = b_n/A_n,n
DOFOR i = n - 1, 1, -1sum = 0
DOFOR col = j + 1, nsum = sum + A_j, col * x_colEND DOx_j = (b_j - sum) / a_j,jEND DOEND Substitute
Error Analysistrue_error = abs(analytical - numerical)cond_number = cond(A,'inf')A_inverse = inv(A)

MATLAB Code

%%% EGR 312%%% Cameron King%%% Coding Assignment 4%%% Problem 1%%% Purpose: The purpose of this code is to solve a system of equations%%% using Gauss elimination with partial pivoting. Then the results will be%%% compared to the backslash operator in an absolute true error%%% calculation. Lastly, the conditioning of the system of equations will%%% be analyzed using the condition number and
%% Set Parameters
A = [4 3 2; 1 3 1; 2 1 3]; % given coefficient matrixb = [960; 510; 610]; % units of copper, zinc, glass respectivelyn = 3; % number of equations
%% Analytical Solution
x0 = A\b; % built-in MATLAB command for solving a system of equations
%% Numerical Solution
% Use Gauss Elimination to simplify to an upper triangular matrixr = 0; % initialize dummy variablefor col = 1:n-1 % variable representing the number of columnsif A(col,col) == 0 % if a diagonal value in the matrix is zerorow = col; % the row number is the same as the column number on the diagonalfor row = row+1:n % variable for the number of rowsendcol1 = A(col,:); % variable for all values in a columnb1 = b(col); % variable for the b value that corresponds with the columnA(col,:) = A(row,:); % reverse the order of the equationsb(col) = b(row); % reverse the order of the equationsA(row,:) = col1; % switch column and rowb(row) = b1; % simplifies the value of bendfor j = 1+r:n-1 % iterative variableC = A(j+1,col)/A(col,col); % variable for the fraction of the leading coefficients that will give the same valueA(j+1,:) = A(j+1,:)-C*A(col,:); % subtract the two equations to get a zerob(j+1) = b(j+1)-C*b(col); % subtract the two b-values as wellendr = r+1; % do this process again for the next elementend
%Back Substitutionx(n) = b(n)/A(n,n); % solves the equation Ax=bfor j = n-1:-1:1 % iterative variable for every equation except the lastsum = 0; % initializing variable in the pseudocodefor col = j+1:n % iterative variable to go through each columnsum = sum+A(j,col)*x(col);endx(j) = (b(j)-sum)/A(j,j); % gives each x valueend
% Display the Solutiondisp('The solution to the system of equations:');disp(x(1)); disp(x(2)); disp(x(3)); % Since x is a row vector, this puts the solutions in column form
%% Error Analysis
A = [4 3 2; 1 3 1; 2 1 3]; % given coefficient matrixformat long % increases the significant digitstrue_error = abs(x0 - x'); % true error between numerical and analyticalformat short % return significant digits back to normal
% Evaluate the condition numbercond_number = cond(A,'inf'); % condition number command
% Check for ill-conditioningA_inverse = inv(A); % inverse of matrix AA_scaled = [1 3/4 1/2; 1/3 1 1/3; 2/3 1/3 1]; % scaled so no number in a row is over 1A_scaledinv = inv(A_scaled); % if values >>1, system is ill-conditioned% largest value is 1.6842 which is not much larger than 1A0 = A_inverse*A; % if not equal to identity matrix, ill-conditioned% A0 equals IA1 = inv(A_inverse); % if not equal to matrix A, ill-conditioned% A1 equals A% Therefore, the system is well-conditioned

Problem 2

Problem Setup

The problem was setup by hand. The Taylor series expansion can be used to derive numerical estimates of higher order derivatives that are called finite difference approximations. Centered different approximations were used in this problem, because they are typically more accurate than forward or backward approximations.

After the system of equations was found, the system was written in matrix notation in MATLAB to find the analytical and numerical solutions.

Analytical Solution

The backslash solver was used to analytical solve this system of equations for the concentration of the chemical. Then, the values for c(0) and c(10) were included to get the entire solution set. Even though the following solution is the analytical solution, there is still a small amount of round-off error due to MATLAB's limited storage space of each value.

The analytical solution, which represents the concentration of the chemical over the distance interval, was plotted to determine trends as distance along the canal is changed. The plot is shown below.

Figure 1: The concentration of the chemical over distance is shown in the plot.

According to the analytical solution and Figure 1, the concentration of the chemical is decreasing as distance is increasing. The slope along each successive step is closer to zero. Given a further distance in the canal, the slope and the concentration may level off at zero.

Numerical Solution

Like in Problem 1, Gauss elimination with partial pivoting was employed to find the numerical solution to the system of equations. The MATLAB results are following.

The MATLAB Code section below shows the exact code that was used to obtain the solution. Notice that the numerical solution is identical to the analytical solution, indicating the effectiveness of Gauss elimination with partial pivoting.

Sample Calculation Using Gauss Elimination With Partial Pivoting

Forward elimination for the first column and the subsequent back substitution is shown in the following sample calculation. In a complete calculation, forward elimination would be completed for each row before using back substitution so that the solutions can be easily found. This sample is intended to show the idea behind Gauss elimination with partial pivoting.

Error Analysis

The error between the analytical and numerical solutions was found using an absolute true error calculation. Since both solutions were the exact same to all the significant digits stored in MATLAB, the absolute true error is a column vector of zeros. The calculation is shown below.

The conditioning of the system of equations was also determined by first calculating the condition number in MATLAB.

cond (A) = 6.0492

Then, the system was checked for ill-conditioning using the same methods as in Problem 1. All three methods failed to show ill-conditioning, as indicated in the MATLAB code below. These methods support the small condition number in concluding that the system of equations is well-conditioned. The three methods are outlined below.

Discussion

The given second order differential equation represents a chemical that is diffusing through a one-dimensional canal. Since it results from a steady state mass balance, I would expect the concentration to be decreasing and approaching a constant value. In other words, I would expect the slope to be approaching zero as distance increases. Figure 1 confirms my hypothesis. The results and plot show that the concentration is highest initially and decreases as distance increases.

Just like in Problem 1, Gauss elimination with partial pivoting accurately solved the system of equations to provide the correct concentration along the canal. The absolute true error between the solution found using Gauss elimination with partial pivoting and MATLAB's backslash solver was zero, proving the effectiveness Gauss elimination with partial pivoting. The high accuracy can partially be attributed to the partial pivoting. Partial pivoting avoids division by zero, minimizes round-off error, and improves ill-conditioning if present. Gauss elimination without partial pivoting should not be employed to solve a complex system of equations due to these advantages.

The condition number provides the condition or "error" within the system. A condition number indicates the change in the output of a system given a small change in the input. All condition numbers are greater than 1, and the condition number of this system is 6.0492. This small value signifies that the system is well-conditioned. The equations in a well-conditioned system have widely differing slopes so that it is easy to find the point of intersection or solution. Conversely, in an ill-conditioned system, the intersection between equations would be difficult to discover because of similar slopes of the equations. Since this problem is well-conditioned, the solution is likely very accurate.

Pseudocode

Set Parameters D, U, deltax, x, A, b, nD = 2;U = 1;k = 0.2;deltax = 2;x = 0:2:10A = [-4.8 1 0 0; 3 -4.8 1 0; 0 3 -4.8 1; 0 0 3 -4.8]b = [-240; 0; 0; -20]n = 4
Analytical Solutionc0 = A\b;c0 = [80; c0; 20]
Numerical Solutionr=0
DOFOR col = 1, n-1
DOIF A_col,col == 0row = col
DOFOR row = row +1, nENDDOcol1 = A_col,lengthb1 = b_colA_col,length = A_row,lengthb_col = b_rowA_row,length = col1b_row = b1ENDDO
DOFOR j = 1+r,n-1f = A_j+1,col / A_col,colA_j+1,length = A_j+1,length - f * A_col,lengthb_j+1 = b_j+1 - f * b_colENDDOr = r +1ENDDO
Substitutec_n = b_n/A_n,n
DOFOR i = n - 1, 1, -1sum = 0
DOFOR col = j + 1, nsum = sum + A_j, col * c_colEND DOc_j = (b_j - sum) / a_j,jEND DOEND Substitute
Plot ResultsDOPLOT x, concentrationInclude title, ylabel, xlabel
Error Analysistrue_error = abs(c0 - c)cond_number = cond(A,'inf')A_inverse = inv(A)

MATLAB Code

%%% EGR 312%%% Cameron King%%% Coding Assignment 4%%% Problem 1%%% Purpose: The purpose of this code is to solve a system of equations%%% using Gauss elimination with partial pivoting. Then the results will be%%% compared to the backslash operator in an absolute true error formula.%%% Lastly, the conditioning of the system will be analyzed using condition%%% number and matrix inverse properties.
%% Set Parameters
A = [-4.8 1 0 0; 3 -4.8 1 0; 0 3 -4.8 1; 0 0 3 -4.8]; % coefficient matrixb = [-240; 0; 0; -20]; % given b-vectorn = 4; % number of equations
%% Analytical Solution
c0 = A\b; % built-in MATLAB command for solving a system of equationsc0 = [80; c0; 20]; % solution with given values for c(0) and c(10)% Gives the concentration of the chemical
%% Numerical Solution
% Use Gauss Elimination to simplify to an upper triangular matrixr = 0; % initialize dummy variablefor col = 1:n-1 % variable representing the number of columnsif A(col,col) == 0 % if a diagonal value in the matrix is zerorow = col; % the row number is the same as the column number on the diagonalfor row = row+1:n % variable for the number of rowsendcol1 = A(col,:); % variable for all values in a columnb1 = b(col); % variable for the b value that corresponds with the columnA(col,:) = A(row,:); % reverse the order of the equationsb(col) = b(row); % reverse the order of the equationsA(row,:) = col1; % switch column and rowb(row) = b1; % simplifies the value of bendfor j = 1+r:n-1 % iterative variablef = A(j+1,col)/A(col,col); % variable for the fraction of the leading coefficients that will give the same valueA(j+1,:) = A(j+1,:)-f*A(col,:); % subtract the two equations to get a zerob(j+1) = b(j+1)-f*b(col); % subtract the two b-values as wellendr = r+1; % do this process again for the next elementend
%Back Substitutionc(n) = b(n)/A(n,n); % solves the equation Ac=bfor j = n-1:-1:1 % iterative variable for every equation except the lastsum = 0; % initializing variable in the pseudocodefor col = j+1:n % iterative variable to go through each columnsum = sum+A(j,col)*c(col);endc(j) = (b(j)-sum)/A(j,j); % final solution of c valuesend
c = [80; c'; 20]; % Since c is a row vector, this puts the solutions in column form% The given values are also included in the complete solution set
% Display the Solutiondisp('The solution to the system of equations:');disp(c)
%% Error Analysis
A = [-4.8 1 0 0; 3 -4.8 1 0; 0 3 -4.8 1; 0 0 3 -4.8]; % coefficient matrixformat long % increases the significant digitstrue_error = abs(c0 - c); % true error between numerical and analyticalformat short % return significant digits back to normal
% Evaluate the condition numbercond_number = cond(A,'inf'); % condition number command
% Check for ill-conditioningA_inverse = inv(A); % inverse of matrix AA_scaled = [1 -1/4.8 0 0; -3/4.8 1 -1/4.8 0; 0 -3/4.8 1 -1/4.8; 0 0 -3/4.8 1]; % scaled matrix so no number in a row is over 1A_scaledinv = inv(A_scaled); % if values >>1, system is ill-conditioned% largest value is 1.3887 which is not much larger than 1A0 = A_inverse*A; % if not equal to identity matrix, ill-conditioned% A0 equals IA1 = inv(A_inverse); % if not equal to matrix A, ill-conditioned% A1 equals A% Therefore, the system is well-conditioned

Reflection

The peer review process was beneficial for me to understand another perspective on how to solve the problem and receive feedback on how to improve my solution. I first modified my initial submission using the peer review feedback. Initially, I didn't realize that Figure 1 was too large to fit on the screen so I decreased it's size immediately. I also moved the Pseudocode sections for both problems right before the respective MATLAB code sections. Then, I decreased the font size of the pseudocode per my partner's suggestions. Lastly, to make the code and pseudocode easier to read, I included more spaces to separate sections. Overall, the peer review comments were very helpful in pointing out formatting improvements.

From reviewing my partner's code, I learned two main things. First, I learned that code is much easier to read when there is adequate spacing between different sections and functions. This is why I decided to include more spaces in my code. Secondly, I learned that creating the full code within various functions makes it very easy to perform the same function on a new system. I will take this into account when working on future problems in MATLAB, especially when the same method needs to be applied to many systems. If written correctly, it will be faster to solve systems but it also will be more challenging to troubleshoot if it fails.

From the initial submission, I made several changes after looking at my peer review comments and my initial submission rubric. The most dramatic change made was moving the backslash solver to the analytical solution and solving the system using a numerical method in problem 2. I solved it using Gauss elimination with partial pivoting and changed the pseudocode to reflect this. I also included sample calculations for both numerical solutions to show the reader how Gauss elimation with partial pivoting is done by hand. Then, with an analytical and numerical solution, I included an absolute true error calculation to determine the error using Gauss elimination with partial pivoting. I also changed the title and corresponding explanation of the plot from problem 2 to showcase concentration over distance rather than concentration over time. Additionally, I corrected the condition number calculation from problem 1 to use the original coefficient matrix.

Page updated

Report abuse