1Some Computational Science Algorithms and Data StructuresKeshav PingaliUniversity of Texas, AustinComputational science• Simulations of physical phenomena– fluid flow over aircraft (Boeing 777)– fatigue fracture in aircraft bodies– evolution of galaxies–….• Two main approaches– continuous methods: fields and differential equations (eg. Navier-Stokes equations, Maxwell’s equations,…)– discrete methods/n-body methods: particles and forces (eg. gravitational forces)• We will focus on continuous methods in this lecture– most differential equations cannot be solved exactly– must use numerical methods that compute approximate solutions• discretization: convert calculus problem to linear algebra problem– finite-difference, finite-element and spectral methodsOrganization• Finite-difference methods– ordinary and partial differential equations– discretization techniques• explicit methods: Forward-Euler method• implicit methods: Backward-Euler method• Finite-element methods– mesh generation and refinement– weighted residuals• Key algorithms and data structures– matrix computations• algorithms– matrix-vector multiplication (MVM)– matrix-matrix multiplication (MMM)– solution of systems of linear equations» direct methods» iterative methods• data structures– dense matrices– sparse matrices– graph computations• mesh generation and refinementOrdinary differential equations• Consider the odeu‘(t) = -3u(t)+2u(0) = 1• This is called an initial value problem– initial value of u is given– compute how function u evolves for t > 0• Using elementary calculus, we can solve this ode exactlyu(t) = 1/3 (e-3t+2)2/32Problem• For general ode’s, we may not be able to express solution in terms of elementary functions• In most practical situations, we do not need exact solution anyway– enough to compute an approximate solution, provided • we have some idea of how much error was introduced• we can improve the accuracy as needed• General solution: – convert calculus problem into algebra/arithmetic problem• discretization: replace continuous variables with discrete variables• in finite differences, – time will advance in fixed-size steps: t=0,h,2h,3h,…– differential equation is replaced by difference equationForward-Euler method• Intuition:– we can compute the derivative at t=0 from the differential equationu‘(t) = -3u(t)+2– so compute the derivative at t=0 and advance along tangent to t =h to find an approximation to u(h)• Formally, we replace derivative with forward difference to get a difference equation–u’(t) Æ (u(t+h) – u(t))/h• Replacing derivative with difference is essentially the inverse of how derivatives were probably introduced to you in elementary calculusBack to ode• Original odeu‘(t) = -3u(t)+2• After discretization using Forward-Euler:(u(t+h) – u(t))/h = -3u(t)+2• After rearrangement, we get difference equationu(t+h) = (1-3h)u(t)+2h • We can now compute values of u:u(0) = 1u(h) = (1-h)u(2h) = (1-2h+3h2)…..Exact solution of difference equation• In this particular case, we can actually solve difference equation exactly• It is not hard to show that if difference equation isu(t+h) = a*u(t)+bu(0) = 1the solution is u(nh) = an+b*(1-an)/(1-a)• For our difference equation, u(t+h) = (1-3h)u(t)+2hthe exact solution isu(nh) = 1/3( (1-3h)n+2)• Stability:– values computed from difference equation will blow up if • ||(1-3h)|| > 1 Î h > 2/3– for this problem, forward-Euler is stable only if step size is less than 2/3– in general, forward-Euler is stable only for small enough step sizes3Comparison• Exact solutionu(t) = 1/3 (e-3t+2)u(nh) = 1/3(e-3nh+2) (at time-steps)• Forward-Euler solution uf(nh) =1/3( (1-3h)n+2)• Use series expansion to compareu(nh) = 1/3(1-3nh+9/2 n2h2 …+ 2)uf(nh) = 1/3(1-3nh+n(n-1)/2 9h2+…+2)So error = O(nh2) (provided h < 2/3)• Conclusion:– error per time step (local error) = O(h2)– error at time nh = O(nh2)h=1/3h=.2h=0.1h=0.01exact solutionChoosing time step• Time-step needs to be small enough to capture highest frequency phenomenon of interest• Nyquist’s criterion– sampling frequency must be at least twice highest frequency to prevent aliasing– for most finite-difference formulas, you need sampling frequencies (much) higher than the Nyquist criterion• In practice, most functions of interest are not band-limited, so use– insight from application or– reduce time-step repeatedly till changes are not significant• Fixed-size time-step can be inefficient if frequency varies widely over time interval– other methods like finite-elements permit variable time-steps as we will see latertimeBackward-Euler method• Replace derivative with backward differenceu’(t+h) Æ (u(t+h) – u(t))/h• For our ode, we getu(t+h)-u(t)/h = -3u(t+h)+2which after rearrangementu(t+h)= (2h+u(t))/(1+3h)• As before, this equation is simple enough that we can write down the exact solution:u(nh) = ((1/(1+3h))n+ 2)/3 • Using series expansion, we getu(nh) = (1-3nh + (-n(-n-1)/2) 9h2+ ...+2)/3u(nh) = (1 -3nh + 9/2 n2h2+ 9/2 nh2+...+2)/3So error = O(nh2) (for any value of h)h=1000h=0.1h=0.01exact solutionComparison• Exact solutionu(t) = 1/3 (e-3t+2)u(nh) = 1/3(e-3nh+2) (at time-steps)• Forward-Euler solution uf(nh) =1/3( (1-3h)n+2)error = O(nh2) (provided h < 2/3)• Backward-Euler solution ub(n*h) = 1/3 ((1/(1+3h))n+ 2)error = O(nh2) (h can be any value you want)• Many other discretizationschemes have been studied in the literature– Runge-Kutta– Crank-Nicolson– Upwind differencing–…Red: exact solutionBlue: Backward-Euler solution (h=0.1)Green: Forward-Euler solution (h=0.1)4Systems of ode’s• Consider a system of coupled ode’s of the formu'(t) = a11*u(t) + a12*v(t) + a13*w(t) + c1(t)v'(t) = a21*u(t) + a22*v(t) + a23*w(t) + c2(t)w'(t) = a31*u(t) + a32*v(t) + a33*w(t) + c3(t)• If we use Forward-Euler method to discretizethis system, we get the following system of simultaneous equationsu(t+h)–u(t) /h = a11*u(t) + a12*v(t) + a13*w(t) + c1(t)v(t+h)–v(t) /h = a21*u(t) + a22*v(t) + a23*w(t) + c2(t)w(t+h)–w(t) /h= a31*u(t) + a32*v(t) + a33*w(t) + c3(t)Forward-Euler (contd.)• Rearranging, we getu(t+h) = (1+ha11)*u(t) + ha12*v(t) + ha13*w(t) + hc1(t)v(t+h) = ha21*u(t) + (1+ha22)*v(t) + ha23*w(t) + hc2(t)w(t+h) = ha31*u(t) + ha32*v(t) + (1+a33)*w(t) + hc3(t)• Introduce vector/matrix notationu(t) = [u(t)
View Full Document
Unlocking...