Engineering Computation Laboratory: Solving Ordinary Differential Equations (ODEs) - Prof., Lab Reports of Materials science

Download Engineering Computation Laboratory: Solving Ordinary Differential Equations (ODEs) - Prof. and more Lab Reports Materials science in PDF only on Docsity! 10/24/2007 1 MAT E 460 ENGINEERING COMPUTATION LABORATORY INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS Prof. Antonios Zavaliangos LeBow 441, x2078, azavalia@coe.drexel.edu 10/24/2007 2 Where do you find ODEs? MOST PROBLEMS IN ENGINEERING INVOLVE ODE 10/24/2007 5 Numerical Solutions Few differential equations can be solved exactly. Numerical solutions are desirable. Key questions: - What is the level of accuracy desired? - What is the computational cost? We will focus on step-by-step methods (also called difference or discrete variable methods) 10/24/2007 6 Euler’s method ii ii xx yy dx dy − − ≈ + + 1 1It is the simplest - it uses the fact that ( ) ),(),( 11 yxfxxyyyxfdx dy iiii −+=⇔= ++ Two possibilities: Euler Forward: (explicit) Euler Backward (implicit): ( ) ),(11 iiiiii yxfxxyy −+= ++ ( ) ),( 1111 ++++ −+= iiiiii yxfxxyy SIMPLE MAY LEAD TO NON-LINEAR EQUATION 10/24/2007 7 Forward vs. Backward Backward Forward Known starting point 10/24/2007 10 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 dx=1. dx=0.5 dx=0.2 dx=0.1 Analytical Example: Euler Backward (implicit) 1)0(,5 =−≈ yy dx dy Euler backward is UNCONDITIONALLY STABLE (which compensates for complexity) 10/24/2007 11 Forward or Backward? Backward Forward If the accuracy here is OK then forward is better Backward is preferred when it provides adequate accuracy for increments below the stability limit of the forward More steps Error 10/24/2007 12 Runge-Kutta method(s) Attempts to improve accuracy over Euler’s methods by evaluating the derivative at subintervals of each step ( ) 21 121 2 , 2 ;, kxyy kxyxxfkyxfk ii iiii ⋅Δ+= ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ Δ + Δ +== + RK order 4 ( ) ( ) ( ) ( ) ( ) 6/22 ,;2/,2/ 2/,2/;, 43211 3423 121 kkkkxyy xkyxxfkxkyxxfk xkyxxfkyxfk ii iiii iiii +++Δ+= Δ+Δ+=Δ+Δ+= Δ+Δ+== + RK order 2 ),( yxf dx dy = 10/24/2007 15 How to solve a system of linear ODES There are many methods - extensions of those used for a single equation. e.g., Euler forward is shown below: xttttxxtyy ttxygyy ttxyfxx ttxtyg tt yy ttxtyf tt xx ttxtyg dt tdy ttxtyf dt tdx nnnnnn nnnnn nnnnn nnn nn nn nnn nn nn ,),(),(where, ),,( ),,( )),(),(( )),(),(( )),(),(()( )),(),(()( 1 1 1 1 1 1 1 −=Δ== Δ+= Δ+= ⇔ = − − = − − ⇔ = = + + + + + + + ttxygyy ttxyfxx nnnnn nnnnn Δ+= Δ+= ++++ ++++ ),,( ),,( 1111 1111 Euler backward shown below leads to 2 non-linear algebraic equations with 2 unknowns Extensions to other methods (R-K, Predictor-Corrector etc.) are possible 11, ++ nn xy 10/24/2007 16 How to solve a system of non-linear algebraic equations n nn n nnnn n nn n nnnn y fyy x fxxyxfyxf y fyy x fxxyxfyxf yxf yxf ∂ ∂ −+ ∂ ∂ −+=≈ ∂ ∂ −+ ∂ ∂ −+=≈ ⇔ = = −− −− 2 1 2 122 1 1 1 111 2 1 )()(),(0),( )()(),(0),( 0),( 0),( Shown below is the solution for a 2 non-linear system of algebraic equations. Generalization to n equations is possible. ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ − ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ − + + ),( ),( 2 1 1 22 11 1 1 nn nn nn nn n n n n yxf yxf y f x f y f x f y x y x The matrix to be inverted is called the Jacobian 10/24/2007 17 Application: Molecular Dynamics • Computes the motion of individual particles in models of solids, liquids, and gases. • Simulation is based on the integrations of differential equations for motion of N-atoms • Contrary to Monte-Carlo simulation it is a deterministic method (pseudo-stochastic aspects can be incorporated - Langevin dynamics (thermal reservoir)) •See tutorial at http://www.mse.ncsu.edu/CompMatSci/Tutorial/PS/ov37.html and simulation example at http://www.ims.uconn.edu/centers/simul/progdoc/atomdemo/atomdemo.html 10/24/2007 20 How to solve a system of linear ODES There are many methods - extensions of those used for a single equation. e.g., Euler forward is shown below: xttttxxtyy ttxygyy ttxyfxx ttxtyg tt yy ttxtyf tt xx ttxtyg dt tdy ttxtyf dt tdx nnnnnn nnnnn nnnnn nnn nn nn nnn nn nn ,),(),(where, ),,( ),,( )),(),(( )),(),(( )),(),(()( )),(),(()( 1 1 1 1 1 1 1 −=Δ== Δ+= Δ+= ⇔ = − − = − − ⇔ = = + + + + + + + ttxygyy ttxyfxx nnnnn nnnnn Δ+= Δ+= ++++ ++++ ),,( ),,( 1111 1111 Euler backward shown below leads to 2 non-linear algebraic equations with 2 unknowns Extensions to other methods (R-K, Predictor-Corrector etc.) are possible 11, ++ nn xy 10/24/2007 21 How to solve a system of non-linear algebraic equations n nn n nnnn n nn n nnnn y fyy x fxxyxfyxf y fyy x fxxyxfyxf yxf yxf ∂ ∂ −+ ∂ ∂ −+=≈ ∂ ∂ −+ ∂ ∂ −+=≈ ⇔ = = −− −− 2 1 2 122 1 1 1 111 2 1 )()(),(0),( )()(),(0),( 0),( 0),( Shown below is the solution for a 2 non-linear system of algebraic equations. Generalization to n equations is possible. ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ − ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ = ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ − + + ),( ),( 2 1 1 22 11 1 1 nn nn nn nn n n n n yxf yxf y f x f y f x f y x y x The matrix to be inverted is called the Jacobian 10/24/2007 22 Molecular dynamics ∑= ))(),(),(),((1)( 2 2 txtxtxtxF mdt tdx jijiij i && )()( ))(),(),(),((1)( tv dt tdx tvtvtxtxF mdt tdv i i jijiij i = = ∑ Is equivalent to: If N atoms are considered we have 6N equations (3 for position and 3 for velocity for 3D simulations 10/24/2007 25 BOTTOM LINE • Solving ODEs is key to many engineering problems • MANY methods exist • Selection should be based on - accuracy - stability - computation cost 10/24/2007 26 EXTRA INFO The rest 3 slides are from material that I will not cover this time. You do not need to know it – it is here just for reference. 10/24/2007 27 Other options : Multi-step methods Improve the accuracy of the single step methods by using information from more than one prior points ),(where( simply of instead ,, 1 1 1 1 11 iii iii kk m k kik m k kiki yxff xfyy tscoefficienbafbxyay = Δ+= =⋅Δ+= + = −+ = −++ ∑∑