Steady-State Mass Balance

Prompt

Linear algebraic equations arise in the solution of differential equations. For example, the following 2^ndorder homogeneous differential equation results from a steady-state mass balance for a chemical in a one-dimensional canal (see figure 1 below), where 𝑐 = concentration, 𝑡 = time, 𝑥 = distance, 𝐷 = diffusion coefficient, 𝑈 = fluid velocity, and𝑘 = a first-order decay rate. Convert this differential equation to an equivalent system of simultaneous algebraic equations. Given 𝐷 = 2, 𝑈 = 1, 𝑘 = 0.2, 𝑐(0) = 80, and 𝑐(10) = 20, solve these equations from 𝑥 = 0 to 10 with Δ𝑥 = 2, and develop a plot of concentration versus distance. What would you expect to see? What do your results show? Show all of your work, especially for how you set up the problem as a linear system. (Hint: Refer to Chapter 4 where we discussed finite-divided-difference approximations to derivatives. Determine a general approximation for the equation and then develop a system of equations using the Δ𝑥 = 2 step size from 𝑥 = 0 to 10.).

Figure 1: Equation used for concentration

Problem Setup (for both approaches)

We can use finite-divided-difference approximations of the derivatives to solve for c. The objective is to convert the provided equation into a system of equations to solve for concentration at various distances. The solution steps and example calculation are provided simultaneously below. Once the system of equations is constructed, MATLAB can be used to easily solve for the concentration vector 'c' via Naive Gauss (numerical approach) or the backslash solver (analytical approach).

Numerical Solution Steps

Below we use Naive Gauss elimination to solve the system. This method involves eliminating variables in equations 2 through n (where n is total number of rows) by subtracting preceding rows that are scaled by some number. We do this until each row n has its pivot column equal to n. The last variable x(n) can be solved for directly in row n and can then be substituted back into preceding equations to solve for the remaining variables.

See below for solution steps with Naive Gauss using a 3x3 coefficient matrix (the solution steps are identical for our four by four setup).

Pseudocode for Numerical Approach

%% Variable declarations and input data

% Declare matrix A (coefficients) and B (of form Ax = b)

a = ...;

b = ...;

% Collect size information about A. 'n' is used for number of equations.

[n,m] = size(a);

% Allocate space for solution vector 'x'

x = size(n);

%% Forward elimination loops

    DO for k = 1:n - 1

        for i = k + 1: n

            % Find value to multiply row n by

            factor = a(i,k) / a(k,k);

            % Subtract the altered row n from row n + 1 by navigating each column individually.

            % Now the leading variable in row n+1 has a coefficient of zero

            DO for j = k + 1: n

                a(i,j) = a(i,j) - factor*a(k,j);

            end DO

            % Subtract scaled value of b(n) in row n from b(n+1) in row n+1.

            b(i) = b(i) - factor*b(k);

        end DO

    end DO

%% Back substitution operations

    % Solve for last variable since there is only one unknown.

    x(n) =  b(n)/a(n,n);

    % Begin substitution loop. Starts at second to last row and works its

    % way up the matrix

    DO for i = n - 1: -1: 1

        % Temporary variable sum declared

        sum = b(i);

        % Use sum to isolate variable being solved for in each row

        DO for j = i + 1:n

            sum = sum - a(i,j)*x(j);

        end DO

        % Divide a variable's corresponding coefficient over to get the

        % variable alone on one side. This determines the variable's value

        x(i) = sum/a(i,i);

    end DO

Analytical Solution Steps

The analytical solution is conducted using MATLAB's backslash solver. Inputting the following matrix into the command window or a script file to use the backslash function returns the solution vector in the workspace.

Pseudocode for Analytical Approach

% Declare x vector

% Declare coefficient matrix

% Declare B matrix of form Ax = b

% Solve for concentration using backslash solver

c = A\B;

% Add outer bounds c(0) and c(10)

c = [80 c' 20];

% Plot function

Figure 2: Concentration vs distance plot

Results

Figure 3: Results for concentration from both approaches

The numerical and analytical approaches resulted in identical solution vectors, as seen in figure 3.

Plotting the resulting solution (identical between methods) with respect to distance shows that the chemical's concentration decreases with increasing distance, and begins to level off the closer it gets to zero. This is how we would expect a chemical to behave when dispered in a fluid. Naturally chemicals move from high concentrations to low concentrations as they diffuse into the environment. As distance increases, however, the concentration begins to approach zero at a slower rate. This is to be expected because at significant distances there will still theoretically be some concentration of the chemical (although this isn't always true in real-life applications).

Error Analysis

Looking at the condition number of the coefficient matrix A can give us a sense of error. Ill-conditioned systems refer to systems of equations where small changes to the coefficient matrix lead to significant changes in the solution matrix x (or in this case, 'C'). These systems may introduce error, especially in the form of round-off error and truncation error, because a large set of potential solutions to the system will satisfy the equations. As we perform row operations, scale, and eliminate variables from the matrix, any amount of rounding (which the computer must do if there are enough decimals on a coefficient) may lead to inaccurate solutions. You can think of an ill-conditioned system as one whose solution matrix is very sensitive to changes in A.

Condition numbers are used for coefficient matrices to determine whether they are ill-conditioned or well-conditioned. Having a high condition number can mean your system is ill-conditioned and is prone to truncation and round-off errors. Having a condition number close to 1 means the system is well-conditioned and we can be confident about the solution to a greater number of significant figures.

Our matrix A has a condition number of 4.7309, calculated using MATLAB's built-in function cond(A) (see appendix I), meaning our system is well-conditioned. We can calculate the number of lost significant figures using the following equation:

= log(condition #)

= log(4.7309)

= 0.675

This small rounding error originates from our low condition number and supports our solution for c. Just because our system of equations is well-conditioned does not necessarily mean our solutions are "accurate". Because our numerical and analytical approaches yeilded the same results, we can expect to see no error between them. See the figure below for the comparison between the two approaches.

Figure 4: Percent error between analytical and numerical approaches

Appendix I: Numerical approach: MATLAB code

clear all; close all; clc; %% Variable declarations and input data % Declare matrix A (coefficients) and B (of form Ax = b)a = [-1.2,0.25,0,0;0.75,-1.2,0.25,0;0,0.75,-1.2,0.25;0,0,0.75,-1.2];b = [-60;0;0;-5]; % Collect size information about A. 'n' is used for number of eqns.[n,m] = size(a); % Allocate space for solution vector 'x'x = size(n); %% Forward elimination loops for k = 1:n - 1 for i = k + 1: n % Find value to multiply row n by factor = a(i,k) / a(k,k); % In coefficient matrix. Subtracts the altered row n % from row n + 1 by navigating each column individually. % Now the leading variable in row n+1 has a coefficient of zero for j = k + 1: n a(i,j) = a(i,j) - factor*a(k,j); end % In matrix 'b'. Subtracts scaled value of b(n) in row n from % b(n+1) in row n+1. This preserves validity of system (must % operate on 'both sides') b(i) = b(i) - factor*b(k); end end %% Back substitution operations % Solves for last variable (and in this case in the last row) directly % since there is only one unknown. x(n) = b(n)/a(n,n); % Begin substitution loop. Starts at second to last row and works its % way up the matrix for i = n - 1: -1: 1 % Temporary variable used to hold and transport the known aspects % of each row (each equation) to the same side as 'b'. Initialize % this variable with the b(i) value of its row sum = b(i); % Uses sum to move the individual products of all known variables and their % coefficients over to the side of 'b' in each row. This isolates the % unknown variable and prepares it to be evaluated by dividing both % sides by its coefficient for j = i + 1:n sum = sum - a(i,j)*x(j); end % Divides a variable's corresponding coefficient over to get the % variable alone on one side. This determines the variable's value x(i) = sum/a(i,i); end % Add bounds. Display x = [80 x' 20]; display(x);

Appendix II: Analytical approach: MATLAB code

clear all; close all; % Declare x vectorx = 0:2:10; % Coefficient matrix & display condition numberA = [-1.2,0.25,0,0;0.75,-1.2,0.25,0;0,0.75,-1.2,0.25;0,0,0.75,-1.2];display(cond(A)); %Solution matrix of form Ax = bB = [-60;0;0;-5]; % Solve for concentration using backslash solverc = A\B; % Add outer bounds c(0) and c(10)c = [80 c' 20]; % Plot functionfigure(1)hold all;plot(x,c,'-o', 'LineWidth',3,'Markersize',2,'Color',[17 17 17]/255);xline(0); yline(0);xlabel('Distance');ylabel('Concentration');title('Concentration with Increasing Distance');xlim([0 10]);legend('c(x)','location', 'northeast')box on; grid on;

Google Sites

Report abuse