Scientific Computing

1. Scientific Computing Numerical Solution Of Ordinary Differential Equations - Further Analysis

2. Solving First Order ODE’s • A general form for a first order ODE is • Or alternatively • We desire a solution x(t) which satisfies this ODE and one specified boundary condition. x(a) = c

3. Numerical Solution • The value of the true solution is approximated at a set of n values of t. • We denote the approximation at these pts by so that

4. Numerical Solution • Let x(t) be the actual solution for dx/dt at the step values. • Then, assuming no roundoff error, the difference in the calculated and true value is the truncation error,

5. Taylor Series • Last time we looked at expressing the solution x(t) about some starting point t0 using a Taylor expansion. • The ODE is given by

6. Error Analysis for Taylor Approximation • Then, if we use the first n terms of the series as an approximation for the solution, we get:

7. Error Analysis for Taylor Approximation • We saw that the local truncation error is bounded as follows: • where

8. Truncation Error vs Roundoff Error

9. Stability of solutions Solution of ODE is • Stable if solutions resulting from perturbations of initial value remain close to original solution • Asymptotically stable if solutions resulting from perturbations converge back to original solution • Unstable if solutions resulting from perturbations diverge away from original solution without bound

10. Example: Stable Solutions

11. Example: Unstable solution • ODE x’ = x is unstable • (solution is x(t)=cet ) • we show solutions with Euler’s method

12. Example: Unstable solution

13. Example: asymptotically stable solutions • ODE x’ = -x is stable • (solution is x(t)= ce-t ) (solution -> 0) • if h too large, numerical solution is unstable • we show solutions with Euler’s method in red

14. x’=-x, stable but slow solution

15. x’=-x, stable, rather inaccurate solution for h large

16. x’=-x, stable but poor solution

17. Stability Analysis • The solution to satisfying some initial condition x(a) = s will be: • unstable if fx ≥ δ for some positive δ and • stable if fx ≤ -δ for some positive δ. • Proof: Let x(t,s) be the solution to the IVP above. Then Since t and s are independent variables, then dt/ds=0, so ->

18. Stability Analysis • Proof(continued): So, if we assume x is C1 Let The solution to this diff eq is u(t) = c eQ(t) , where If q(r) = fx(r,x(r,s)) ≥ δ >0 then Q-> ∞ so, u-> ∞ and the solution is unstable. If q(r) = fx(r,x(r,s)) ≤ -δ <0 then Q-> -∞ so, u-> 0 and the solution is stable. QED

19. Euler’s Method • Review: • The simple Euler method is a linear approximation, and is a perfect solution only if the function is linear (or at least linear in the interval). • This is inherently inaccurate. • Because of this inaccuracy small step sizes are required when using the algorithm.

20. Better Methods • If we took more terms in the Taylor’s Series expansion of f(t,x), we would get a more accurate solution algorithm. • Problem: It is difficult to compute higher-order derivatives numerically. • Idea: Approximate the value of higher derivatives of f(t,x) by evaluating f several times between iterates: ti and ti+1 .

21. Better Methods • The ODE +IV x(a)= x0 =s can be alternatively solved as an integral equation. For example, on the first interval from x0 to x1 :

22. Better Methods In general: • The problem is that we don’t know x(t), so evaluating the integral is not trivial. • We assume that if h is small enough, f(t,x) will not change that much over the interval [ti, ti+1]. Using “left-rectangular integration,” • This is the Euler Method. A family of improved versions are called Runge-Kutta methods

23. Euler Method Euler Solution x xi true solution h ti ti+1 = ti + h Tuncation error: true solution

24. Euler vs Actual Solution Euler Solution RK methods differ in how they estimate fi For Euler Method (a type of RK method), fi ≈ f(ti, xi) x xi true solution fi h ti ti+1 = ti + h

25. Runge-Kutta Order 2 Method x xi fi h Weighted average of two slopes x ti ti+αh ti+1