Introduction to Computer Science • Sedgewick and Wayne • Copyright © 2007 • http://www.cs.Princeton.EDU/IntroCS9. Scientific ComputingApplications of Scientific ComputingScience and engineering challenges.Fluid dynamics.Seismic surveys.Plasma dynamics.Ocean circulation.Electronics design.Pharmaceutical design.Human genome project.Vehicle crash simulation.Global climate simulation.Nuclear weapons simulation.Molecular dynamics simulation.Common features.Problems tend to be continuous instead of discrete.Algorithms must scale to handle huge problems.Commercial applications.Web search.Financial modeling.Computer graphics. Digital audio and video.Natural language processing.Architecture walk-throughs.Medical diagnostics (MRI, CAT).Floating PointIEEE 754 representation.Used by all modern computers.Scientific notation, but in binary.Single precision: float = 32 bits.Double precision: double = 64 bits.Ex. Single precision representation of -0.453125.1 0 1 1 1 1 1 0 1 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0sign bitexponent significand125 1/2 + 1/4 + 1/16 = 0.8125-1-1 × 2125 - 127 × 1.8125 = -0.453125biasphantom bitFloating PointRemark. Most real numbers are not representable, including π and 1/10.Roundoff error. When result of calculation is not representable. Consequence. Non-intuitive behavior for uninitiated.Financial computing. Calculate 9% sales tax on a 50¢ phone call.Banker's rounding. Round to nearest integer, to even integer if tie.if (0.1 + 0.2 == 0.3) { // NO }if (0.1 + 0.3 == 0.4) { // YES }double a1 = 1.14 * 75; // 85.49999999999999 double a2 = Math.round(a1); // 85double b1 = 1.09 * 50; // 54.50000000000001 double b2 = Math.round(b1); // 55SEC violation (!)you lost 1¢Catastrophic CancellationA simple function.Goal. Plot f(x) for -4 ⋅ 10-8 ≤ x ≤ 4 ⋅ 10-8.Exact answer€ f (x) = 1 − cos xx2Catastrophic CancellationA simple function.Goal. Plot f(x) for -4 ⋅ 10-8 ≤ x ≤ 4 ⋅ 10-8.€ f (x) = 1 − cos xx2IEEE 754 double precision answerCatastrophic CancellationEx. Evaluate fl(x) for x = 1.1e-8.Math.cos(x) = 0.99999999999999988897769753748434595763683319091796875.(1.0 - Math.cos(x)) = 1.1102e-16(1.0 - Math.cos(x)) / (x * x) = 0.9175Catastrophic cancellation. Devastating loss of precision when small numbers are computed from large numbers, which themselves are subject to roundoff error.nearest floating point value agrees withexact answer to 16 decimal places.inaccurate estimate of exact answer (6.05 ⋅ 10-17)80% larger than exact answer (about 0.5)public static double fl(double x) { return (1.0 - Math.cos(x)) / (x * x);}Numerical CatastrophesAriane 5 rocket. [June 4, 1996]10 year, $7 billion ESA project exploded after launch.64-bit float converted to 16 bit signed int.Unanticipated overflow.Vancouver stock exchange. [November, 1983]Index undervalued by 44%.Recalculated index after each trade by adding change in price.22 months of accumulated truncation error.Patriot missile accident. [February 25, 1991]Failed to track scud; hit Army barracks, killed 28.Inaccuracy in measuring time in 1/20 of a secondsince using 24 bit binary floating point.Copyright, ArianespaceGaussian Elimination0 x0+1 x1+ 1 x2= 4 2 x0+4 x1- 2 x2= 2 0 x0+3 x1+ 15 x2= 36 € A =0 1 12 4 −20 3 15⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ , b =4236⎡ ⎣ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ Linear System of EquationsLinear system of equations. N linear equations in N unknowns.Fundamental problems in science and engineering.Chemical equilibrium.Linear and nonlinear optimization.Kirchoff's current and voltage laws.Hooke's law for finite element methods.Leontief's model of economic equilibrium.Numerical solutions to differential equations.…matrix notation: find x such that Ax = bChemical EquilibriumEx. Combustion of propane.Stoichiometric constraints.Carbon: 3x0 = x2.Hydrogen: 8x0 = 2x3.Oxygen: 2x1 = 2x2 + x3.Normalize: x0 = 1.Remark. Stoichiometric coefficients tend to be small integers;among first hints suggesting the atomic nature of matter.x0C3H8 + x1O2 ⇒ x2CO2 + x3H2OC3H8 + 5O2 ⇒ 3CO2 + 4H2Oconservation of massKirchoff's Current LawEx. Find current flowing in each branch of a circuit.Kirchoff's current law.10 = 1x0 + 25(x0 - x1) + 50 (x0 - x2).0 = 25(x1 - x0) + 30x1 + 1(x1 - x2).0 = 50(x2 - x0) + 1(x2 - x1) + 55x2.Solution. x0 = 0.2449, x1 = 0.1114, x2 = 0.1166.x0x1x2conservation of electrical chargeUpper Triangular System of EquationsUpper triangular system. aij = 0 for i > j.Back substitution. Solve by examining equations in reverse order.Equation 2: x2 = 24/12 = 2.Equation 1: x1 = 4 - x2 = 2.Equation 0: x0 = (2 - 4x1 + 2x2) / 2 = -1.for (int i = N-1; i >= 0; i--) { double sum = 0.0; for (int j = i+1; j < N; j++) sum += A[i][j] * x[j]; x[i] = (b[i] - sum) / A[i][i];} 2 x0+ 4 x1- 2 x2= 2 0 x0+ 1 x1+ 1 x2= 4 0 x0+ 0 x1+ 12 x2= 24Gaussian EliminationGaussian elimination.Among oldest and most widely used solutions.Repeatedly apply row operations to make system upper triangular.Solve upper triangular system by back substitution.Elementary row operations.Exchange row p and row q.Add a multiple α of row p to row q.Key invariant. Row operations preserve solutions.pq0 x0+1 x1+ 1 x2= 4 2 x0+4 x1- 2 x2= 2 0 x0+3 x1+ 15 x2= 36 2 x0+4 x1- 2 x2= 2 0 x0+1 x1+ 1 x2= 4 0 x0+3 x1+ 15 x2= 36 2 x0+ 4 x1- 2 x2= 2 0 x0+ 1 x1+ 1 x2= 4 0 x0+ 0 x1+ 12 x2= 24 (interchange row 0 and 1)(subtract 3x row 1 from row 2)Gaussian Elimination: Row OperationsElementary row operations.Forward elimination. Apply row operations to make upper triangular.Pivot. Zero out entries below pivot app.for (int p = 0; p < N; p++) { for (int i = p + 1; i < N; i++) { double alpha = A[i][p] / A[p][p]; b[i] -= alpha * b[p]; for (int j = p; j < N; j++) A[i][j] -= alpha * A[p][j]; }}€ aij= aij−aipappapjbi= bi−aipappbpGaussian Elimination: Forward EliminationppForward elimination. Apply row operations to make upper triangular.Pivot. Zero out entries below pivot app.Gaussian Elimination: Forward Elimination€ * * * * ** * * * ** * * * ** * * * ** * * * *⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥