CSE 245 ComputerAided Circuit Simulation and Verification Lecture

CSE 245: Computer-Aided Circuit Simulation and Verification Lecture Note 5. 1 Numerical Integration Prof. Chung-Kuan Cheng 1

Numerical Integration: Outline • • Richardson extrapolation (Bulirsch-Stoer) Rosenbrock method (Runge Kutta) Predictor-corrector method Matrix exponential 2

Integration Methods 3

Bulirsch-Stoer Method Caveats: • Nonsmooth function: RK • Contain singular points: RK • Very smooth and right-hand sides expensive to compute: Predictor-corrector 4

Bulirsch-Stoer Method Approach: • Modified midpoint method • Extrapolation • Stepsize control 5

Bulirsch-Stoer Method: midpoint method Given dx/dt=f(t, x), H and n, set h=H/n z 0=x(t) z 1=z 0+hf(t, z 0) zm+1=zm-1+2 hf(t+mh, zm) for m=1, 2, . . , n-1 x(t+H): xn=1/2[zn+zn-1+hf(t+H, zn)] Error: xn-x(t+H)=∑i=1 aih 2 i [1, 2, 3] Example Sequence: n=2, 4, 6, 8, 10, 12, 14, …(Deuflhard) 6

Bulirsch-Stoer Method: Extrapolation T 00 T 10 T 20 … T 11 T 21 … T 22 … … Tk 0=xk Tk, j+1=Tkj+(Tkj-Tk-1, j)/[(nk/nk-j)2 -1], j=0, 1, …, k-1 Solution: Tkk Error: |Tkk-Tk, k-1| Errk: H 2 k+1 7

Bulirsch-Stoer Method: Stepsize Control Stepzie Hk=HS 1(S 2/errk)1/(2 k+1) Complexity A 0=n 0+1 Ak+1=Ak+nk+1 Work per unit step Wk=Ak/Hk Strategy minimize Wk ([4]. 17. 3. 3) For y(x+H)≈yn+(yn-yn/2)/3, we use 1. 5 derivative evaluations per step h. For Runge-Kutta, it takes 4 evaluations. 8

Rosenbrock Methods ODE: dx/dt=f(x) Stepsize: h Process: x(t 0+h)=x 0+∑i=1, sbiki (1 -rhf’)ki=hf(x 0+∑j=1, i-1 aijkj)+hf’∑j=1, i-1 rijkj, i=1, …, s Runge-Kutta: r=rij=0 for all ij. 9

Rung-Kutta Method (4 th order) For 4 th order RK method, we evaluate the derivatives four times: once at the initial points, twice at trial midpoints, and once at a trial endpoint. The final solution is calculated from the 4 derivatives. k 1= hf(tn, xn) k 2=hf(tn+0. 5 h, xn+0. 5 k 1) k 3= hf(tn+0. 5 h, xn+0. 5 k 2) k 4= hf(tn+h, xn+k 3) xn+1=xn+1/6 k 1+1/3 k 2+1/3 k 3+1/6 k 4+O(h 5) 10

Predictor-Corrector Methods: Admas-Bashforth ODE: dx/dt=f(x) Predictor xn+1=xn+h/12(23 fn-1 -16 fn-1+5 fn-2)+O(h 4) Corrector xn+1=xn+h/12(5 fn+1+8 fn-fn-1)+O(h 4) 11

References 1. J. A. Gaunt, The deferred approach to the limit, IIinterpenetrating lattices, Trans. Roy, Soc. , Lond. 226, 350 -361, 1927 2. R. De Vogelaere, On a paper of Gaunt concerned with the start of numerical solutions of differential equations, Z. Angew. Math. Phys, 151 -156, 1957 3. W. B. Gragg, On extrapolation algorithms for ordinary initial value problems, J. of SIAM, 384 -403, 1965 4. W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes, 3 rd Edition, 2007 12