Scientific computing
This presentation is the property of its rightful owner.
Sponsored Links
1 / 36

Scientific Computing PowerPoint PPT Presentation


  • 75 Views
  • Uploaded on
  • Presentation posted in: General

Scientific Computing. Linear Systems – Gaussian Elimination. Linear Systems. Linear Systems. Solve Ax=b, where A is an n  n matrix and b is an n  1 column vector We can also talk about non-square systems where A is m  n , b is m  1, and x is n  1

Download Presentation

Scientific Computing

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.


- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript


Scientific computing

Scientific Computing

Linear Systems – Gaussian Elimination


Linear systems

Linear Systems


Linear systems1

Linear Systems

  • Solve Ax=b, where A is an nnmatrix andb is an n1 column vector

  • We can also talk about non-square systems whereA is mn, b is m1, and x is n1

    • Overdetermined if m>n:“more equations than unknowns”

    • Underdetermined if n>m:“more unknowns than equations”


Singular systems

Singular Systems

  • A is singular if some row islinear combination of other rows

  • Singular systems can be underdetermined:or inconsistent:


Using matrix inverses

Using Matrix Inverses

  • To solve Ax = b it would be nice to use the inverse to A, that is, A-1

  • However, it is usually not a good idea to compute x=A-1b

    • Inefficient

    • Prone to roundoff error


Gaussian elimination

Gaussian Elimination

  • Fundamental operations:

    • Replace one equation with linear combinationof other equations

    • Interchange two equations

    • Multiply one equation by a scalar

  • These are called elementary row operations.

  • Do these operations again and again to reduce the system to a “trivial” system


Triangular form

Triangular Form

  • Two special forms of matrices are especially nice for solving Ax=b:

  • In both cases, successive substitution leads to a solution


Triangular form1

Triangular Form

  • A is lower triangular


Triangular form2

Triangular Form

  • Solve by forward substitution:


Triangular form3

Triangular Form

  • Solve by forward substitution:


Triangular form4

Triangular Form

  • Solve by forward substitution:

Etc


Triangular form5

Triangular Form

  • If A is upper triangular, solve by back-substitution:


Triangular form6

Triangular Form

  • Solve by back-substitution:

Etc


Gaussian elimination algorithm

Gaussian Elimination Algorithm

  • Do elementary row operations on the augmented system [A|b] to reduce the system to upper triangular form.

  • Then, use back-substitution to find the answer.


Gaussian elimination1

Gaussian Elimination

  • Example:

  • Augmented Matrix form:


Gaussian elimination2

Gaussian Elimination

  • Row Ops: What do we do to zero out first column under first pivot?

  • Zero out below second pivot:


Gaussian elimination3

Gaussian Elimination

  • Back-substitute


Gaussian elimination4

Gaussian Elimination


Matlab implementation

Matlab Implementation

  • Task: Implement Gaussian Elimination (without pivoting) in a Matlab M-file.

  • Notes

    • Input = Coefficient matrix A, rhs b

    • Output = solution vector x


Matlab implementation1

Matlab Implementation

  • Class Exercise: We will go through this code line by line, using the example in Pav section 3.3.2 to see how the code works.


Matlab implementation2

Matlab Implementation

  • Matlab:

    • >> A=[2 1 1 3; 4 4 0 7; 6 5 4 17; 2 -1 0 7]

    • A =

    • 2 1 1 3

    • 4 4 0 7

    • 6 5 4 17

    • 2 -1 0 7

    • >> b = [7 11 31 15]'

    • b =

    • 7

    • 11

    • 31

    • 15

>> gauss_no_pivot(A,b)

ans =

1.5417

-1.4167

0.8333

1.5000


Matlab implementation3

Matlab Implementation

  • Class Exercise: How would we change the Matlab function gauss_no_pivot so we could see the result of each step of the row reduction?


Gaussian elimination pivoting

Gaussian Elimination - Pivoting

  • Consider this system:

  • We immediately run into a problem:we cannot zero out below pivot, or back-substitute!

  • More subtle problem:


Gaussian elimination pivoting1

Gaussian Elimination - Pivoting

  • Conclusion: small diagonal elements are bad!

  • Remedy: swap the row with the small diagonal element with a row below, this is called pivoting


Gaussian elimination pivoting2

Gaussian Elimination - Pivoting

  • Our Example:

  • Swap rows 1 and 2:

  • Now continue:


Gaussian elimination pivoting3

Gaussian Elimination - Pivoting

  • Two strategies for pivoting:

    • Partial Pivoting

    • Scaled Partial Pivoting


Partial pivoting

Partial Pivoting


Matlab partial pivoting

Matlab – Partial Pivoting

  • Partial Pivoting:

    • At step k, we are working on kth row, pivot = Akk . Search for largest Aik in kth column below (and including) Akk . Let p = index of row containing largest entry.

    • If p ≠ k, swap rows p and k.

    • Continue with Gaussian Elimination.


Matlab partial pivoting1

Matlab – Partial Pivoting

  • Finding largest entry in a column:

    • >> A

    • A =

    • 2 1 1 3

    • 4 4 0 7

    • 6 5 4 17

    • 2 -1 0 7

    • >> [r,m] = max(A(2:4,2))

    • r =

    • 5

    • m =

    • 2

  • Why isn’t m = 3?


Matlab partial pivoting2

Matlab – Partial Pivoting

  • Swapping rows m and k:

    • BB =

    • 2 1 1 3

    • 4 4 0 7

    • 6 5 4 17

    • 2 -1 0 7

    • >> BB([1 3],:) = BB([3 1], :)

    • BB =

    • 6 5 4 17

    • 4 4 0 7

    • 2 1 1 3

    • 2 -1 0 7


Matlab partial pivoting3

Matlab – Partial Pivoting

  • Code change?

  • All that is needed is in main loop (going from rows k ->n-1) we add

    • % Find max value (M) and index (m) of max entry below AAkk

    • [M,m] = max(abs(AA(k:n,k)));

    • m = m + (k-1); % row offset

    • % swap rows, if needed

    • if m ~= k, AA([k m],:) = AA([m k],:);

    • end


Matlab partial pivoting4

Matlab – Partial Pivoting

  • Class Exercise

  • Review example in Moler Section 2.6


Scaled partial pivoting

Scaled Partial Pivoting

Pav, Section 3.2


Practice

Practice

  • Class Exercise: We will work through the example in Pav, Section 3.2.2


Practice1

Practice

  • Class Exercise: Do Pav Ex 3.7 for partial pivoting (“naïve Gaussian Elimination”).

  • Do Pav Ex 3.7 for scaled partial pivoting.


  • Login