Numerical Methods Applied to Chemical Engineering: Sample Homework Solutions

Homework 0 – Problem 3 Commented [U1]: You should title your homework problem. An accurate estimate of the heat capacity (CP) for a particular inlet mixture is critical to an efficient separation, so you decide to carry out a series of heat capacity measurements at different temperatures in the operating range and fit a polynomial model to compute specific heat capacity as a function of temperature. To find the polynomial that best fits the data, we seek the n values of the weights di id 0 that give the best fit (in a least squares sense) for the linear system of nT equations, where nT is the number of T values at which CP(T) was measured: d 0  d1T j  d 2T j 2    d nd T j nd  CP (T j ), j  1,..., nT The measured values of the heat capacity have a relative error bound of 0.03 and absolute error bound of 50 J/(kg-K). Use the 2-norm throughout this problem for consistency. 1. For the case of a square system with nT = nd + 1 there will be values of d that fit the data. Express the system as Ad = c, where d contains the weights. Write out the forms of A and c. Matrix-vector form of the polynomial fit for CP(T) Write the nT equations for specific heat capacity assuming nT = nd + 1: Commented [U2]: Short introductions to the problem statement can be nice to define your nomenclature and set the stage for your answers. However, it is not necessary and if you do this limit your description to be as short as possible. MAKE IT LESS THAN a quarter of a page. Commented [U3]: Adding lines/spacing can alert the grader to the start of a new part of the problem. Again, a nice feature but not necessary. d 0  d1T1  d 2T12    d nd T1nd  CP (T1 ) d 0  d1T2  d 2T2 2    d nd T2 nd  CP (T2 ) d 0  d1T3  d 2T32    d nd T3nd  CP (T3 ) In vector matrix form:  d 0  d1TnT  d 2TnT 2    d nd TnT nd  CP (TnT ) nd  2   d   C (T )  1 T1 T1  T1   0   P 1  nd   2 1 T    T T  2 2 2    d1   CP (T2 )  nd   2 1 T    T3  T3   d 2    CP (T3 )  3                   n  2 1 TnT TnT  TnT d   d nd  CP (TnT ) From this, the forms of A and c can easily be determined: 1 T1 T12  T1nd   CP (T1 )      n 2 1 T  CP (T2 )  T2  T2 d  2    c   CP (T3 )  A  1 T3 T32  T3nd                C (T ) n 2 1 TnT TnT  TnT d   P nT  1 Commented [U4]: Boxing your final answers and giving short introductory sentences can be very helpful for graders. I suggest you do this. 2. Let the T measurements be evenly spaced between 273 K and 373 K, using nd = 2,3,…,8 (with nT = 3,4,…,9) plot cond(A) and norm(A-1) on the same axes as a function of nT. Discussion of the algorithm to produce A, cond(A), and norm(A-1) The paulson_HW2_P2.m file submitted online accesses a user-written function (condA_norminvA_errorPredCp()), reference for more information regarding the specifics of the algorithm) that generates the matrix A for each specified size (3 through 9). Using A, the function returns both its condition number and the norm of its inverse using the built-in MATLAB®functions cond() and norm(), respectively. The condA_norminvA_errorPredCp() function generates the matrix A using two nested for loops, the first of which spans the rows and the second of which spans the columns. The first row (i = 1) corresponds to the lower bound for temperature (T_low). The second loop creates entries in matrix A (1 to nT) that correspond to the proper power of temperature (i.e. T0 to TnT-1). After each row is filled with the correct elements of T along its columns, the temperature value (T) is adjusted to an evenly spaced value between T_low and T_high (through the addition step of T_spacing). This step occurs before the row value, i, is updated, which is important because each row in A must correspond to a different temperature value. The algorithm repeats this process until increment i reaches the final row in A (i = nT). Since the problem statement asks for parameters at various nT values, I decided to add a third nested loop above the two that generate A which runs from a lower to an upper bound on nT. After the two loops (discussed above) generate A, the condition number, norm of A inverse, and absolute error of CP (see part 4) are stored in vectors based on the corresponding nT value. The outer loop repeats this process from a user-input lower bound of nT (nT_low) to a user-input upper bound of nT (nT_high) After the three loops are completely iterated, these parameter vectors (that have stored every specified parameter for all nT’s input to the function) are designated as the function output. These vectors (cond_A, norm_invA, abs_error_Cp_pred) are then plotted versus nT in the main function. Condition number of A and norm of A inverse The condition number of A and the norm of A inverse are plotted versus nT below on Figure P2.1. The algorithm used to generate this plot is described above (reference paulson_HW2_P2.m for more information regarding any part of problem 2). Figure P2.1 shows that the condition number of A grows in order of magnitudes with small changes in nT. Additionally, their linear trends on a semilog plot imply that these parameters have an exponential correlation with respect to nT. 2 Commented [U5]: Here is an example of going TOO FAR. I describe the submitted code in detail. This is not necessary especially for this problem where we want you to worry more about the equations than any algorithm (the algorithm is very simple, just make an A matrix and compute the values). It is completely fine to do this, but a waste of time. I would recommend given a much shorter description (e.g., “My submitted MATLAB code generates matrix A from part 1 and computes the condition number of norm of A inverse using built-in functions”). Commented [U6]: The plot is commented well and nicely labeled. You should do this in your reports. Notice that the quality of the plot is not very good, you can achieve much better quality in a variety of ways (please ask Joel for more information on this). Figure P2.1. Condition Number and norm of inverse for A on a log10 scale as a function of nT 3. Again for the square system, compute relative and absolute error bounds on d for several values of nd. Let the T measurements be evenly spaced between 273 and 373 K for the same nT and nd values in part 2. Plot the relative and absolute error bounds as a function of nT. Note: Express your error bounds in terms of the relative and absolute error bounds in heat capacity. Error matrix rearrangement The matrix vector form of the equation for heat capacity is: Commented [U7]: Derivations such as this can be nice and helpful for you. I would recommend that you include it when it is not obvious as you can reference it during exams (you should be allowed to have all your notes and graded HWs). Here, this is fairly obvious so I could have jumped to the boxed equation directly. Ad  c Add a perturbation (Δd) in d that creates a perturbation (Δc) in c A(d  Δd)  c  Δc Distribute A on the LHS Ad  AΔd  c  Δc Notice the original equation Ad = c, substitute this into the equation above c  AΔd  c  Δc AΔd  Δc 3 Induced Norm Inequality Proof The induced norm of matrix A   M N by the vector x   N is:  Ax   A  max   x  NOTE: THROUGHOUT THIS PROBLEM ||A|| implies 2-norm of matrix A By the definition of the maximum Ax x  A Ax  A x , Δx   N : x  0 *Induced norm inequality proof (referenced throughout the following sections) Absolute Error of d In order to calculate the absolute error of d  Δd  , left multiply the above relationship by the inverse of A Δd  A1Δc Take the 2-norm of the above relationship Δd  A1Δc From the induced norm inequality (proved in part 3) Δd  A1 Δc *Absolute error inequality expression Relative Error of d In order to calculate the relative error of d, take the 2-norm of the original matrix equation: Ad  c Divide the relative error inequality (above) by this expression A1 Δc Δd  Ad c Multiply the norm of Ad on both sides of the equation 4 Commented [U8]: This was an example of something given in the notes that you could have referenced. NOT REQUIRED AND A WASTE OF TIME. You would still need to give the final boxed equation but you would not have to go through any details. Δd  A1 Δc c Ad From the induced norm inequality (proved in part 3) Δd  A1 Δc c A d Divide both sides by the norm of d Δd Δc  A 1 A d c *Relative error inequality expression 2-Norm of Δc Δc  nT  Cp 2 i i 1 Expand the summation term Δc  Cp12  Cp2 2    CpnT 2 Since each of the measurement error terms, ΔCPi are equivalent, they can be group together nT times. ΔCP is the absolute measurement error in heat capacity, ΔCP = 50 J/kg-K (Given). Δc  nT Cp 2 Δc  nT  Cp J kgK *2-Norm of the heat capacity error vector Δc  nT  50 Error bounds for d The given absolute and relative error bounds for the measured heat capacities (c) are: Absolute Error: Δc  nT  50 J kgK Relative Error: 5 Δc c  0.03 The absolute and relative errors for weights d can be calculated using the error inequalities derived above in part 3: Δd Δc  A 1 A Δd  A1 Δc d c Δd  nT  50 J A 1 kgK Δd  0.03 A1 A d Define the norm of A times the norm of its inverse to be the condition number (cond(A)). Δd  nT  50 J A 1 kgK Δd  0.03  cond ( A) d *Used to calculate absolute and relative error bounds of d in terms of the error bounds of c Plot of the error bounds of d as a function of nT The relative and absolute errors of CP weights (d) are plotted as a function of nT below on Figure P2.2. The methodology used to generate these error values for d is described above. (Figure P2.2 also shows cond(A) and the norm(A-1) from part 2 for comparison) Figure P2.2. Relative and absolute error of d plotted on a log10 scale as a function of nT Figure P2.2 shows that the relative error bound of d follows the same trend as the condition number of A and the absolute error bound of d follows the same trend as the 2-norm of A inverse. Due to these trends the error bounds of d grow exponentially fast making this polynomial fit unsuitable for CP as a function of T. 6 4. Again for the square system, derive an expression bounding the absolute uncertainty of the specific heat capacity prediction from the fitted d for some arbitrary 𝑇𝑇� between 273 K and 373 K, for arbitrary nd, and justify the bound. Plot the numerical values for the bounds for the values of nd and nT above, for 𝑇𝑇�= 300 K. Assume: 1 T^   d0     d1      d 2   CPPr ed (T^ )  CPTrue(T^ )  CP(T^ )         d   nd  ^ ^ T2  T nd And that there exists some 𝐝𝐝̂ within the error bound on d computed above such that: 1 T^  ^ T2 ^  T nd  d^0  ^ d  1   d^2       ^  d nd       C True(T^ ) P      Subtracting these two equations yields: 1 T^  ^ T2 ^  T nd   d 0   d^     0   d   d^  1   1    d 2    d^2                 ^   d nd   d nd           C (T^ ) P          �. Define Δd to be the difference between d and 𝐝𝐝̂ and, T to be the polynomial expansion of T Substitute these relations into the above equation ^ Td  CP(T ) Take the 2-norm of both sides of the equation ^ Td  CP(T ) Using the matrix induced norm inequality: Td  T d 7 ^ CP(T )  T d Absolute error inequality from part 3: d  A1 Δc ^ CP(T )  T A1 Δc Since ΔCP is a scalar its 2-norm is equivalent to the absolute value of the quantity itself ^ CP(T )  T A1 Δc Where ΔCP is the difference in the predicted and true heat capacity values � values substituted into the polynomial fit for CP T is the row vector of T A is the temperature matrix representing the CP polynomial fit (part 1) Δc is the absolute measurement error bound for heat capacity = 50 J/(kg-K) Δc  nT  Δc (reference part 3 for derivation) *Used to calculate absolute uncertainty in predicted CP values � = 300 K Predicted CP uncertainty at 𝐓𝐓 The absolute uncertainty in the predicted CP value (defined as |ΔCP|) plotted versus nT = 3,..,9 is shown below in Figure P2.3. The values of this plot were obtained using the condA_norminvA_errorPredCp() function in the submitted paulson_HW2_P2.m file. The algorithm is discussed in detail in part 2. Figure P2.3. Absolute uncertainty in predicted CP at 300 K (log10 scale) 8 Commented [U9]: Here is a nice example of a completely summarized equation. I am telling the graders everything they need to know for the final answer (including defining my terms) while referencing above for derivations. If you define things previously, it is ok to NOT define them again, but remember graders are human and make mistakes. � = 300 K, the predicted specific heat capacity has an absolute Figure P2.3 shows that at T error bound of order 30 which corresponds to a polynomial fit with nine temperature measurements (nT = 9). These large magnitudes in the predicted CP values are a direct result of the large-normed T vector and the ill-conditioned A matrix. Together, these parameters significantly amplify the small absolute error bound in the measured heat capacities. 5. For the more general case nT > nd + 1, derive expressions for absolute error bounds for the leastsquares estimates dLS using the relation seen in Homework 0: AT Ad LS = AT c Derivation Starting with the general least-squares estimate equation: AT Ad LS = AT c Add a perturbation (ΔdLS) in dLS that creates a perturbation (Δc) in c AT A(d LS  Δd LS )  AT (c  Δc) Distribute ATA on the LHS and AT on the RHS AT Ad LS  AT AΔd LS  AT c  AT Δc Notice the original equation AT Ad LS = AT c , cancel on both sides of the equation AT AΔd LS  AT Δc *Least squares (LS) error equation Multiply both sides of the LS error equation by the inverse of ATA: Δd LS  ( AT A)1 AT Δc Take the 2-norm of both sides of this equation Δd LS  ( AT A)1 AT Δc Using similar ideas to the induced norm inequality (proved in part 4), the RHS norm can be shown to be less than the product of the individual vectors/matrices norms. This results in: Δd LS  ( AT A)1 AT Δc *Notes: (1) ||ΔdLS|| is the relative error bound for the least square estimates (dLS) (2) A and c are equivalent to those defined in part 1 9 Commented [U10]: I would recommend giving little descriptions of the plots such as this one so that the plot is not freefloating in space. You can keep it short and sweet. 6. Explain the intuitive meaning of the ill-conditioning of this system For a linear equation, Ax = b, the condition number of matrix A is said defined as the operator of two norms, cond(A) = ||A|| ||A-1||. The condition number is important when solving linear equations because it is an indicator as to how much a perturbation in b will alter the solution vector x. If the condition number is high, small changes in b can create large errors in x. Such a system is said to be ill-conditioned. This usually implies A has small-valued eigenvalues making it close to singular (zero-valued determinant) and/or the system is badly scaled. In this problem, the A matrix equals (part 1): 1 T1  1 T 2  A  1 T3     1 TnT T12 T2 2 T32  TnT 2  T1nd   T2 nd   T3nd       TnTnd  This matrix has a large magnitude discrepancy within each of its columns due to the “Taylorseries type” polynomial fit for heat capacity (CP) in terms of temperature (T) which has increasing orders of T from 0 to nd. Since the values of temperature are usually order 3, these polynomial factors result in extremely large magnitude differences from column to column. This intuitively means that either the chosen scale for the problem parameters are wrong or the fit for CP as a function of temperature was a poor choice. For this particular problem, I believe the polynomial expansion fit was a poor choice since high order functions usually result in large magnitude terms (Terms >= 1010). Additionally, these values span so many orders that it is difficult to estimate each term accurately. Moreover, the large-magnitude terms result in an enormous condition number for A which, in turn, takes the relatively small errors in c and amplifies them within the solution vector d. I would recommend choosing a better function to fit CP with respect to T instead of forcing data into random polynomial expansions and preforming least squares analysis blindly. This more in-depth analysis would require investigating how these variables normally interact with each other in various relationships. Hopefully, this newer fit would eliminate the large magnitude discrepancy between columns; however, rescaling the variables to some small finite range could also be a viable option (i.e. try using relationships that make the temperature and heat capacity dimensionless, maybe through the addition of other important parameters). 7. Use Chebyshev polynomials in place of the Tis to construct matrix A. You can use the MATLAB function ChebyshevPoly. For better conditioning, scale the temperature so that the Chebyshev polynomial is evaluated on [-1, 1] by evaluating at (2T - Tlow - Thigh) / (Thigh - Tlow) rather than at T. Comment on the difference between the error bounds for Tis vs. Chebyshev polynomials on the rescaled temperature range. 10 Discussion of the Chebyshev polynomial algorithm The paulson_HW2_P2.m file submitted online accesses a user-written function (chebyshev_A()), reference for more information regarding the specifics of the algorithm) that generates a modified matrix A for each specified size (3 through 9). Using the ChebyshevPoly() MATLAB function a polynomial is generated for every power of T in the assumed fit for CP. This polynomial is then evaluated at a modified T value on the scale of [-1, 1]. The main idea behind this type of modification is to alter A so that its condition number is significantly lower effectively reducing the error bounds for the solution vector. The userwritten chebyshev_A() function is just a slightly modified version of the cond_A_norminvA_errorPredCp() function discussed in detail in part 2 (reference for more information regarding the framework of this algorithm). The only change to this is inside the third nested for loop, which no longer stores some power of T to the matrix A, but instead performs this substitution/evaluation for the chebyshev polynomial discussed in the paragraph above. From this, it can be seen that a multiple polynomials (evaluated along the altered T scale) are stored in A which effectively adds more entries to the spanning set of the matrix A. Additionally, the rescaled temperature reduces the extremely large magnitude values dealt with in part 2 (i.e. ~2739). Plot of the error bounds of d as a function of nT evaluated using Chebyshev polynomials The relative and absolute errors of CP weights (d) calculated using Chebyshev polynomials are plotted as a function of nT below on Figure P2.4. Figure P2.4. Absolute and relative error bounds for d calculated using Chebyshev polynomials 11 Error bound comparison for d For easier comparison, the relative and absolute errors of d calculated using the full polynomial expansion (part 2) and Chebyshev polynomials (part 7) are plotted together on Figure P2.5 below. The full polynomial fit produces an absolute error range for d of ~104 to ~1010 (over nT = 3 to 9) whereas the Chebyshev polynomials produce a range of ~101.75 to ~102.5. Furthermore, the full polynomial fit produces a relative error range for d of ~105 to ~1026 (over nT = 3 to 9) whereas the Chebyshev polynomials produce a range of ~10-1.5 to ~10-0.5. This shows that the both the absolute and relative error bounds for d drop significantly when Chebyshev polynomials are used in place of the full polynomial expansion to calculate A. Figure 2.5. Error bounds for d using the full polynomial expansion and Chebyshev polynomials 12