Edwards-Penney Chapter 2, sections 2.4, 2.5, 2.6 Topics, Definitions, Theorems, Examples 2.4 Euler's numerical method for y'=f(x,y) ==== ALGORITHM. Euler's Method Given y'=f(x,y), y(x0)=y0 and target x* = x0+(n+1)h, step size h. Find data points (xk,yk) as approximations to solution curve data points (xk,y(xk)) as follows: 1. x[0]=x0, y[0]=y0 [from the initial condition y(x0)=y0] x[1]=x[0]+h, y[1]=y[0]+hf(x[0],y[0]) x[2]=x[1]+h, y[2]=y[1]+hf(x[1],y[1]) x[3]=x[2]+h, y[3]=y[2]+hf(x[2],y[2]) ... x[n+1]=x[n]+h, y[n+1]=y[n]+hf(x[n],y[n]) EXAMPLE. Find the exact solution to y'=x+y/5, y(0)=-3. Then find y(5). ANSWER: y=22 exp(x/5)-5x-25 by the linear integrating factor method. Then y(5)=9.8022002. EXAMPLE. Apply Euler's method to y'=x+y/5, y(0)=-3 with target x* = 5, step size h=1. ANSWERS: Pairs (0,-3), (1,-3.6), (2,-3.32), (3,-1.984), (4,0.6912), (5,4.7430). EXAMPLE. Apply Euler's method to y'=x+y/5, y(0)=-3 with target x* = 1, step size h=0.2. ANSWERS: Pairs (0,-3), (0.2,-3.12), (0.4,-3.205), (0.6,-3.253), (0.8,-3.253), (1,-3.234). EXAMPLE. Falling baseball. Given v'=32-0.16v, v(0)=0, find data points with target x* = 10, step size h=1. ANSWER: (0,0), (1,32), (2,59), (3,81), (4,100), (5,116), (6,130), (7,141), (8,150), (9,158), (10,165). DEF. LOCAL ERROR is the difference between the exact solution to the initial value problem y'=f(x,y), y(x[n])=y[n] evaluated at x=x[n+1] and the Euler approximation y[n+1]=x[n]+hf(x[n],y[n]). DEF. CUMULATIVE ERROR at step n+1 is the difference between the exact solution to the initial value problem y'=f(x,y), y(x[0])=y[0] evaluated at x=x[n+1] and the Euler approximation y[n+1]=x[n]+hf(x[n],y[n]). DEF. ROUNDOFF ERROR is a computer-generated error due to finitely many digits of accuracy. Computer algebra systems can increase the number of digits used beyond hardware limitations, but still there are a finite number of digits available in numerical work. ERROR PROPAGATION is studied in numerical analysis. EXAMPLE. Solve the logistic initial value problem y'=(8-y)y/3, y(0)=1. ANSWER: y=8/(1+7 exp(-8x/3) EXAMPLE. Sensitivity to step size in Euler's method. The logistic model y'=(8-y)y/3, y(0)=1 requires step size adjustments for reasonable results. EXAMPLE. Oscillatory solutions and Euler's method. The problem y'=(cos(x))y, y(0)=1 has periodic solution y=exp(sin(x)), poorly approximated by Euler's method. EXAMPLE. Approximations near vertical asymptotes. The problem y'=x^2+y^2, y(0)=1 is approximated well near x=0 but shows instability near x=1. Computer numerical studies can overlap a singular point in the solution (where y=infinity) thereby giving answers where no answer exists. MAPLE TUTOR for NUMERICAL METHODS # y'=-2xy, y(0)=2, by Euler, Heun, RK4 with(Student[NumericalAnalysis]): InitialValueProblemTutor(diff(y(x),x)=-2*x*y(x),y(0)=2,x=0.5); # The tutor compares exact and numerical slutions. 2.5 Convergence of Euler's method, Heun's method (improved Euler) ==== THEOREM 1. Under smoothness assumptions for the solution y(x) of the problem y'=f(x,y), y(x0)=y0, the Euler appoximations for step size h satisfy an inequality |y[n]-y(x[n])| <= C h where C depends only on f, x0, y0, and the interval [a,b] on which the solution y(x) is defined. SUMMARY: The cumulative error in Euler's method is of order h. EXAMPLE. Apply Euler's method to y'=2xy/(1+x^2), y(0)=1 to estimate the value of C in the theorem. METHOD. Start with h=0.04 and target x* =1.0. Compute C=|y(1)-y[n+1]|/h where x[n+1]=0+nh=x*. Half the step size and repeat. The answer is C is about 0.11. ALGORITHM. The improved Euler method or Heun's method Given y'=f(x,y), y(x0)=y0 with target x*=x0+(n+1)h and step size h. Improved Euler or Heun's method computes approximate data points along the solution curve y(x) by the formulas x[0]=x0, y[0]=y0 (from initial condition y(x0)=y0) x[1]=x[0]+h, y[1]=y[0]+h(f(x[0],y[0])+f(x[1],y*))/2 where y* = y[0]+hf(x[0],y[0]) x[2]=x[1]+h, y[2]=y[1]+h(f(x[1],y[1])+f(x[2],y*))/2 where y* = y[1]+hf(x[1],y[1]) x[3]=x[2]+h, y[3]=y[2]+h(f(x[2],y[2])+f(x[3],y*))/2 where y* = y[2]+hf(x[2],y[2]) ... x[n+1]=x[n]+h, y[n+1]=y[n]+h(f(x[n],y[n])+f(x[n+1],y*))/2 where y* = y[n]+hf(x[n],y[n]) AVERAGE. The factor (f(x[n],y[n])+f(x[n+1],y*))/2 is an average that originates with the integral approximation method known as the Trapezoidal rule. THEOREM. The cumulative error of the Improved Euler method is of order h^2, that is, |y(x[n])-y[n]| <= C h^2. EXAMPLE. Solve y'=x+y, y(0)=1 and evaluate y(1). ANSWER: y=2 exp(x) - x - 1, y(1)=3.4365637. EXAMPLE. Apply Improved Euler to y'=x+y, y(0)=1 with target x* = 1.0 and step size h=0.1. Compare to the Euler method for step sizes h=0.1 and h=0.005. ANSWER: See the table Figure 2.5.4, section 2.5. The table shows Improved Euler getting an extra digit of accuracy for the same step size. Euler's method catches up with step size h=0.005. EXAMPLE. Solve y'=-2xy/(1+x^2), y(0)=1 and evaluate y(1). Answer: y=1/(1+x^2), y(1)=1/2. EXAMPLE. Find the Improved Euler approximations for the problem y'=-2xy/(1+x^2), y(0)=1, with target x* = 1.0 and step size h=0.04. Repeat for smaller step sizes and estimate the error constant C in the theorem |y(x[n])-y[n]| <= C h^2. ANSWER: See Figure 2.5.6, section 2.5. The answer is C=0.12. EXAMPLE. Apply Improved Euler to y'=(8-y)y/3, y(0)=1 with target x* = 5.0 and various step sizes. ANSWER: See Figure 2.5.7 and 2.5.8, section 2.5. The results show poor agreement for coarse steps and excellent agreement with 200 subintervals (h=5/200). EXAMPLE. Apply Improved Euler to y'=(cos(x))y, y(0)=1 with target value x* = 6 Pi and 50, 100, 200 steps. Compare with the exact periodic solution y=exp(sin(x)). ANSWER: See 2.5.9, section 2.5. Excellent results obtained for 100, 200. EXAMPLE. Periodic harvesting logistic problem y'=y(1-y)- sin(2 Pi t) Plot approximate solution curves on a direction field using Improved Euler to generate the threaded curves. ANSWER: See Figure 2.5.12, section 2.5. 2.6 Runge-Kutta numerical method for y'=f(x,y) ==== ALGORITHM. The Runge-Kutta 4th order formulas We assume problem y'=f(x,y), y(x0)=y0 with target x* = x0+(n+1)h and step size h. x[0]=x0, y[0]=y0 from initial condition y(x0)=y0 x[n+1]=x0+nh, k1=f(x[n],y[n]) k2=f(x[n]+h/2,y[n]+h k1/2) k3=f(x[n]+h/2,y[n]+h k2/2) k4=f(x[n]+h,y[n]+h k3) y[n+1]=y[n]+h(k1 + 2 k1 + 2 k2 + k4)/6 EXPLANATION: Derived from Simpson's rule on int(f(x,y(x)),x) in a historical manner, proof uses Taylor series and is not intuitive. The terms are products of predictor-corrector ideas applied to the terms arising from Simpson's rule. THEOREM. The cumulative error of the Runge-Kutta method is of order 4, which means the error satisfies |y(x[n])-y[n]| <= C h^4. EXAMPLE. Apply RK4 to y'=x+y, y(0)=1 with target x* = 1.0 using step size h=0.5. Compare with the exact solution y=2 exp(x) - x - 1 at x=1. ANSWER: The exact value is y(1)=3.4365637. The data pairs from h=0.5 are (0,1), (0.5,1.7969), (1.0,3.4347). For the data pairs with step size h=0.1, see Figure 2.6.1, section 2.5. See Figure 2.6.2 for a comparison with the Improved Euler method. EXAMPLE. Apply RK4 to y'=x^2+y^2, y(0)=1 with x* near the singular point 0.969811. Then compute y(0.5) and its RK4 value with various step sizes. ANSWER: Figure 2.6.4 has the RK4 results, which show first digit errors at x=0.9 when h=0.1. Figure 2.6.5 shows the RK4 results for smaller step sizes, getting consistent results for x=0.9. EXAMPLE. Apply RK4 to y'=x^2+y^2, y(0.9)=14.3049 with x* near the singular point 0.969811. Then compute y(0.95) and its RK4 values with various step sizes. ANSWER: Figure 2.6.6 gives the RK4 results. Agreement is 3 digits minimum with h=0.002, 0.001, 0.0005. The exact solution from Maple is y(0.95) = 50.47230931 using the code MAPLE ANSWER CHECK de:=diff(y(x),x)=x^2+y(x)^2; ic:=y(0.9)=14.3049; dsolve({de,ic},y(x)); subs(x=0.95,%); evalf(%); EXAMPLE. Skydiver problem v' = 9.8 - 0.00016(100 v + 10 v^2 + v^3), v(0)=0. Find the terminal speed 35.578 m/s by roots of equations. Find the speed by RK4 methods and display the results in a table for t from 0 to 20 seconds. ANSWER: Set v'=0 ans solve for v=35.578 by computer. For RK4, use step sizes h=0.2 and h=0.1 looking for consistent results. Figure 2.6.8 shows RK4 results for times 0 to 20 seconds. EXAMPLE. Stiff differential equation. Given y'=5y-6exp(-x), y(0)=1 with exact positive solution y=exp(-x), apply RK4 with h=0.2, 0.1, 0.05 to compute an approximate value for y at target x* = 4.0. ANSWER: Values for the three step sizes at x* = 4.0 are -103, -8903, -631 resp. Since y(4.0)>0, all are wrong with a large error. DEF. An equation y'=f(x,y) is called STIFF if numerical methods applied to it are unstable, giving wrong results or widely varying accuracy depending on step size. Intuition comes from chemical equations where stiffness corresponds to co-existence of very slow and very fast chemical reactions. REMARK. Numerical analysis experts admit that no precise definition of stiffness is possible. HISTORY. The term has been attributed to the mechanical problem x'' + 1001 x' + 1000 x = 0, x(0)=1, x'(0)=0 which represents a damped non-oscillatory spring-mass system. The exact solution is x(t) = -0.001001001001 exp(-1000 t)+1.001001001 exp(-t) It is a stiff spring problem (Hooke's constant 1000) with numerical instability due to heavy damping (1001), seen in sensitivity to step size. MAPLE STIFF EQUATION EXPERIMENT de:=diff(x(t),t,t)+1001*diff(x(t),t)+1000*x(t)=0; ic:=x(0)=1,D(x)(0)=0; q:=dsolve({de,ic},x(t));X:=unapply(rhs(q),t); [X(0),X(0.1),X(0.2)]; # Exact answers p:=dsolve({de,ic},numeric,method=classical[rk4],x(t), stepsize=0.01,output=array([0,0.1,0.2]); ======= end