2nd question in the MatLab section
Assignment 1 Page 1 of 4 Assignment 3 (Linear Systems) Due: Wednesday, March 27, 2023 Theoretical Problems: (due at the beginning of class on the above date) a) The solutions to all theoretical problems are to be neatly handwritten or printed. b) Present problems in the same order as that of the assignment sheet. c) Show all steps for all problems. d) Make sure to start each new problem number at the top of a page. e) STAPLE all problems for HW3 Theoretical Part 1 into one packet with this title at the top of the first page; STAPLE all problems for HW3 Theoretical Part 2 into a separate packet with this title at the top of the first page. (You will turn in two stapled packets.) Unstapled papers will NOT be accepted. f) Also, make sure that your name is on the first page of each packet. MATLAB Problems: (due by 11:45pm on the above date) a) All MATLAB programs will be completed in MATLAB Grader. b) Instructions for each problem are also given in MATLAB Grader. c) Use double precision (default in MATLAB Grader) in all your computer calculations unless otherwise indicated. d) The core of your programs are to be documented (use “%” before comments) Theoretical Problems Part 1: 1) For the MATLAB code, 1: 1: : 2 ( , ) ( , ) ( , )* ( , ) / ( , ) 1: 1 ( , ) ( , ) ( , )* ( , ) for j n for i j for k j n A j j A j j A j k A j k A j j end end for i j A i j A i j A i k A j k end end = = = − = − = + = + Compute the exact number of flops (arithmetic operations) in the entire algorithm in terms of n. Write your final answer as a simplified polynomial, from highest to lowest power. Assume all operations result in defined answers. Page 2 of 4 2) Find the LUPA = factorization (using partial pivoting) of the matrix 1 8 6 1 4 5 2 4 6 A − − − = − − − . Write out P, L, and U for your final answer. 3) Let 1 8 6 1 4 5 2 4 6 A − − − = − − − . Verify that properties (i) and (ii) of the matrix ∞ -norm (i.e., ∞ ⋅ ) are satisfied for A. Now, define another matrix, 3 0 5 1 7 2 0 3 4 B − = − − . Verify that property (iii) is satisfied for both A and B. Theoretical Problems Part 2: 4) Consider the following linear system 1 2 2 1 1 2 1 1 x x ε ε− − − = − , where the computed solution is + = 1 1 ε cx and 10 < ε . use the 1-norm (i.e., 1 ⋅ ) for the following: a) find the condition number in terms of ε . b) evaluate the condition number from part (a) at 510ε −= ; write your answer in scientific notation rounded to one significant digit. using this result, at most how many digits could be trusted in any computed solution, in general (not just cx ), if the exact solution were not known? c) now, using the computed solution cx above, find the relative backward error in terms of ε . d) without finding the exact solution, use your results from parts (a) and (c) to find an upper bound for the relative (forward) error in the solution cx in terms of ε . e) evaluate the upper bound from part (d) at 510ε −= ; write your answer in scientific notation rounded to one significant digit. f) now, using the exact solution and cx , find the true relative forward error in terms of ε . (you can compute the exact solution by hand.) g) evaluate the true relative forward error from part (f ) at 510ε −= ; write your answer in scientific notation rounding to one significant digit. compare its size with that of your answer from part (e). h) compute [ ] 0 lim cond( )a ε→ where 2 1 2 1 a ε− − = − . what happens to the matrix a when 0=ε ? 5) using your results from matlab problem 2 part (c): a) what were the backward errors for ge and gepp? (scientific notation rounded to 1 significant digit is sufficient.) b) what were the forward errors for ge and gepp? (scientific notation rounded to 1 significant digit is sufficient.) page 3 of 4 c) for ge, were the forward and backward errors very different or similar in size (according to the exponent in the scientific notation)? d) for gepp, were the forward and backward errors very different or similar in size (according to the exponent in the scientific notation)? e) given your answers from parts (c ) and (d), would you say that matrix 2a is ill- or well-conditioned, and why? matlab problems (to be completed in matlab grader): 1) in many applications (e.g., numerical solutions of certain differential equations) the coefficient matrix, a, in an nxn linear system, bxa = , takes a tridiagonal form: 1 1 1 2 2 2 3 3 3 4 4 2 1 1 1 n n n n n v w u v w u v w a u v w u v w u v − − − − = , where only the main diagonal and the diagonals immediately above and below the main diagonal potentially have nonzero elements, and all other elements are 0. in such cases, only data in these 3 diagonals is stored as 3 vectors, rather than an entire matrix: the main diagonal could be stored in an nx1 vector, v, and the other diagonals could be stored in (n-1)x1 vectors, u and w. in the matlab script below, only where indicated, modify the “full matrix” version of gaussian elimination, forward substitution, and backward substitution, so that only the above 3 vectors are used, along with the nx1 vector b, an nx1 vector y, containing the output of forward substitution, and an nx1 vector x, containing the final solution output from backward substitution (6 vectors total). note: • the function, gael(), will perform only gaussian elimination – no row exchanges - and return modified vectors u and v (w will be unchanged); • the forward substitution function, forsub(), will only use u and v to compute y • the backward substitution function, backsub(), will only use v, w, and y to compute x. modify only code within the three above functions, where indicated. 2) in this problem, you will convert a function that performs ge (gaussian elimination) into one that performs gepp (gaussian elimination with partial pivoting). a) at the bottom of the script below are three functions: gael, forsub, and backsub. just above these three functions, add a fourth function, called gaelpp, which is a modification of the function gael. write this new function so that it accepts only the a matrix as input and produces two outputs: one will be the modified (overwritten) version of a, which contains the multipliers in the lower triangular portion and the upper triangular matrix generated by gepp; the other output will be the n × n permutation matrix p also generated during gepp. • you could use the matlab function eye( ) to initialize p. page 4 of 4 • write code to search for the appropriate maximum absolute values. you can use a loop with an if-statement or a matlab function like[m,i] = max(a(j:n,j)), which will return an index i for the location of the maximum value of the column vector a(j:n,j); however, you would need to adjust this index. • to swap rows, it may help to use a variable temp to temporarily store a row vector. • you will also swap rows of p accordingly. • overwrite appropriate entries of a with the multipliers instead of saving them in m. b) above the function you wrote in part (a), where indicated at the top of the script, write a piece of code to test it with 1 3.03 12.1 14 3.03 12.1 7 6.11 14.2 21 a − = − − − and 1 119 120 139 b − = − . • assign a1 and b1. • call the functions forsub and backsub to solve the system for x1 with appropriate inputs (these functions are already given at the bottom of the script; note that a matlab script file can call multiple functions that are defined at the bottom of that script file.) • compute the norm of the residual (backward error), 1 1ax b− , using the matlab function norm( ), which computes the 2-norm. (you do not need to use loops to compute the residual: keep in mind that you can use matrix operators to compute errors. in elearning, see “matlab and matrix operations” on p.9 in the matlab tutorial. also note the backward error from ge already computed below.) c) farther down the script, where indicated, write another piece of code to test your function, gaelpp( ), using the 20x20 matrix 2a where ( )2 , 1 1i j a i j = + − . • you can use nested loops to assign a2 . • let the exact solution to some problem, 2 2a x b= , be the 20x1 vector [111...1] tx = (in matlab, the built-in function ones( ) will be useful) • use matrix multiplication in matlab to compute the right-hand side b2. • now, solve the system for x2 by doing the same thing that you did in part (b). • compute the norm of the residual (backward error), 2 2ax b− . also, compute the norm of the forward error, 2x x− , to compare the known exact solution with the solution you computed via gaelpp(). also, note the backward/forward errors from ge already computed below. ε="" .="" use="" the="" 1-norm="" (i.e.,="" 1="" ⋅="" )="" for="" the="" following:="" a)="" find="" the="" condition="" number="" in="" terms="" of="" ε="" .="" b)="" evaluate="" the="" condition="" number="" from="" part="" (a)="" at="" 510ε="" −=";" write="" your="" answer="" in="" scientific="" notation="" rounded="" to="" one="" significant="" digit.="" using="" this="" result,="" at="" most="" how="" many="" digits="" could="" be="" trusted="" in="" any="" computed="" solution,="" in="" general="" (not="" just="" cx="" ),="" if="" the="" exact="" solution="" were="" not="" known?="" c)="" now,="" using="" the="" computed="" solution="" cx="" above,="" find="" the="" relative="" backward="" error="" in="" terms="" of="" ε="" .="" d)="" without="" finding="" the="" exact="" solution,="" use="" your="" results="" from="" parts="" (a)="" and="" (c)="" to="" find="" an="" upper="" bound="" for="" the="" relative="" (forward)="" error="" in="" the="" solution="" cx="" in="" terms="" of="" ε="" .="" e)="" evaluate="" the="" upper="" bound="" from="" part="" (d)="" at="" 510ε="" −=";" write="" your="" answer="" in="" scientific="" notation="" rounded="" to="" one="" significant="" digit.="" f)="" now,="" using="" the="" exact="" solution="" and="" cx="" ,="" find="" the="" true="" relative="" forward="" error="" in="" terms="" of="" ε="" .="" (you="" can="" compute="" the="" exact="" solution="" by="" hand.)="" g)="" evaluate="" the="" true="" relative="" forward="" error="" from="" part="" (f="" )="" at="" 510ε="" −=";" write="" your="" answer="" in="" scientific="" notation="" rounding="" to="" one="" significant="" digit.="" compare="" its="" size="" with="" that="" of="" your="" answer="" from="" part="" (e).="" h)="" compute="" [="" ]="" 0="" lim="" cond(="" )a="" ε→="" where="" 2="" 1="" 2="" 1="" a="" ε−="" −="" ="" −="" ="" .="" what="" happens="" to="" the="" matrix="" a="" when="" 0="ε" 5)="" using="" your="" results="" from="" matlab="" problem="" 2="" part="" (c):="" a)="" what="" were="" the="" backward="" errors="" for="" ge="" and="" gepp?="" (scientific="" notation="" rounded="" to="" 1="" significant="" digit="" is="" sufficient.)="" b)="" what="" were="" the="" forward="" errors="" for="" ge="" and="" gepp?="" (scientific="" notation="" rounded="" to="" 1="" significant="" digit="" is="" sufficient.)="" page="" 3="" of="" 4="" c)="" for="" ge,="" were="" the="" forward="" and="" backward="" errors="" very="" different="" or="" similar="" in="" size="" (according="" to="" the="" exponent="" in="" the="" scientific="" notation)?="" d)="" for="" gepp,="" were="" the="" forward="" and="" backward="" errors="" very="" different="" or="" similar="" in="" size="" (according="" to="" the="" exponent="" in="" the="" scientific="" notation)?="" e)="" given="" your="" answers="" from="" parts="" (c="" )="" and="" (d),="" would="" you="" say="" that="" matrix="" 2a="" is="" ill-="" or="" well-conditioned,="" and="" why?="" matlab="" problems="" (to="" be="" completed="" in="" matlab="" grader):="" 1)="" in="" many="" applications="" (e.g.,="" numerical="" solutions="" of="" certain="" differential="" equations)="" the="" coefficient="" matrix,="" a,="" in="" an="" nxn="" linear="" system,="" bxa="," takes="" a="" tridiagonal="" form:="" 1="" 1="" 1="" 2="" 2="" 2="" 3="" 3="" 3="" 4="" 4="" 2="" 1="" 1="" 1="" n="" n="" n="" n="" n="" v="" w="" u="" v="" w="" u="" v="" w="" a="" u="" v="" w="" u="" v="" w="" u="" v="" −="" −="" −="" −="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ="" ,="" where="" only="" the="" main="" diagonal="" and="" the="" diagonals="" immediately="" above="" and="" below="" the="" main="" diagonal="" potentially="" have="" nonzero="" elements,="" and="" all="" other="" elements="" are="" 0.="" in="" such="" cases,="" only="" data="" in="" these="" 3="" diagonals="" is="" stored="" as="" 3="" vectors,="" rather="" than="" an="" entire="" matrix:="" the="" main="" diagonal="" could="" be="" stored="" in="" an="" nx1="" vector,="" v,="" and="" the="" other="" diagonals="" could="" be="" stored="" in="" (n-1)x1="" vectors,="" u="" and="" w.="" in="" the="" matlab="" script="" below,="" only="" where="" indicated,="" modify="" the="" “full="" matrix”="" version="" of="" gaussian="" elimination,="" forward="" substitution,="" and="" backward="" substitution,="" so="" that="" only="" the="" above="" 3="" vectors="" are="" used,="" along="" with="" the="" nx1="" vector="" b,="" an="" nx1="" vector="" y,="" containing="" the="" output="" of="" forward="" substitution,="" and="" an="" nx1="" vector="" x,="" containing="" the="" final="" solution="" output="" from="" backward="" substitution="" (6="" vectors="" total).="" note:="" •="" the="" function,="" gael(),="" will="" perform="" only="" gaussian="" elimination="" –="" no="" row="" exchanges="" -="" and="" return="" modified="" vectors="" u="" and="" v="" (w="" will="" be="" unchanged);="" •="" the="" forward="" substitution="" function,="" forsub(),="" will="" only="" use="" u="" and="" v="" to="" compute="" y="" •="" the="" backward="" substitution="" function,="" backsub(),="" will="" only="" use="" v,="" w,="" and="" y="" to="" compute="" x.="" modify="" only="" code="" within="" the="" three="" above="" functions,="" where="" indicated.="" 2)="" in="" this="" problem,="" you="" will="" convert="" a="" function="" that="" performs="" ge="" (gaussian="" elimination)="" into="" one="" that="" performs="" gepp="" (gaussian="" elimination="" with="" partial="" pivoting).="" a)="" at="" the="" bottom="" of="" the="" script="" below="" are="" three="" functions:="" gael,="" forsub,="" and="" backsub.="" just="" above="" these="" three="" functions,="" add="" a="" fourth="" function,="" called="" gaelpp,="" which="" is="" a="" modification="" of="" the="" function="" gael.="" write="" this="" new="" function="" so="" that="" it="" accepts="" only="" the="" a="" matrix="" as="" input="" and="" produces="" two="" outputs:="" one="" will="" be="" the="" modified="" (overwritten)="" version="" of="" a,="" which="" contains="" the="" multipliers="" in="" the="" lower="" triangular="" portion="" and="" the="" upper="" triangular="" matrix="" generated="" by="" gepp;="" the="" other="" output="" will="" be="" the="" n="" ×="" n="" permutation="" matrix="" p="" also="" generated="" during="" gepp.="" •="" you="" could="" use="" the="" matlab="" function="" eye(="" )="" to="" initialize="" p.="" page="" 4="" of="" 4="" •="" write="" code="" to="" search="" for="" the="" appropriate="" maximum="" absolute="" values.="" you="" can="" use="" a="" loop="" with="" an="" if-statement="" or="" a="" matlab="" function="" like[m,i]="max(A(j:n,j))," which="" will="" return="" an="" index="" i="" for="" the="" location="" of="" the="" maximum="" value="" of="" the="" column="" vector="" a(j:n,j);="" however,="" you="" would="" need="" to="" adjust="" this="" index.="" •="" to="" swap="" rows,="" it="" may="" help="" to="" use="" a="" variable="" temp="" to="" temporarily="" store="" a="" row="" vector.="" •="" you="" will="" also="" swap="" rows="" of="" p="" accordingly.="" •="" overwrite="" appropriate="" entries="" of="" a="" with="" the="" multipliers="" instead="" of="" saving="" them="" in="" m.="" b)="" above="" the="" function="" you="" wrote="" in="" part="" (a),="" where="" indicated="" at="" the="" top="" of="" the="" script,="" write="" a="" piece="" of="" code="" to="" test="" it="" with="" 1="" 3.03="" 12.1="" 14="" 3.03="" 12.1="" 7="" 6.11="" 14.2="" 21="" a="" −="" ="" ="" ="−" −="" ="" ="" −="" ="" and="" 1="" 119="" 120="" 139="" b="" −="" ="" ="" ="" ="" ="" −="" ="" .="" •="" assign="" a1="" and="" b1.="" •="" call="" the="" functions="" forsub="" and="" backsub="" to="" solve="" the="" system="" for="" x1="" with="" appropriate="" inputs="" (these="" functions="" are="" already="" given="" at="" the="" bottom="" of="" the="" script;="" note="" that="" a="" matlab="" script="" file="" can="" call="" multiple="" functions="" that="" are="" defined="" at="" the="" bottom="" of="" that="" script="" file.)="" •="" compute="" the="" norm="" of="" the="" residual="" (backward="" error),="" 1="" 1ax="" b−="" ,="" using="" the="" matlab="" function="" norm(="" ),="" which="" computes="" the="" 2-norm.="" (you="" do="" not="" need="" to="" use="" loops="" to="" compute="" the="" residual:="" keep="" in="" mind="" that="" you="" can="" use="" matrix="" operators="" to="" compute="" errors.="" in="" elearning,="" see="" “matlab="" and="" matrix="" operations”="" on="" p.9="" in="" the="" matlab="" tutorial.="" also="" note="" the="" backward="" error="" from="" ge="" already="" computed="" below.)="" c)="" farther="" down="" the="" script,="" where="" indicated,="" write="" another="" piece="" of="" code="" to="" test="" your="" function,="" gaelpp(="" ),="" using="" the="" 20x20="" matrix="" 2a="" where="" (="" )2="" ,="" 1="" 1i="" j="" a="" i="" j="+" −="" .="" •="" you="" can="" use="" nested="" loops="" to="" assign="" a2="" .="" •="" let="" the="" exact="" solution="" to="" some="" problem,="" 2="" 2a="" x="" b="," be="" the="" 20x1="" vector="" [111...1]="" tx="(in" matlab,="" the="" built-in="" function="" ones(="" )="" will="" be="" useful)="" •="" use="" matrix="" multiplication="" in="" matlab="" to="" compute="" the="" right-hand="" side="" b2.="" •="" now,="" solve="" the="" system="" for="" x2="" by="" doing="" the="" same="" thing="" that="" you="" did="" in="" part="" (b).="" •="" compute="" the="" norm="" of="" the="" residual="" (backward="" error),="" 2="" 2ax="" b−="" .="" also,="" compute="" the="" norm="" of="" the="" forward="" error,="" 2x="" x−="" ,="" to="" compare="" the="" known="" exact="" solution="" with="" the="" solution="" you="" computed="" via="" gaelpp().="" also,="" note="" the="" backward/forward="" errors="" from="" ge="" already="" computed="">