MCS 471 Practice Problems 5: Numerical Solution of ODEs

Background Reading:

• Gerald and Wheatley, Chapter 5 and beginning of 6

• 5. Numerical Solution of ODEs
  • 5.3 Euler Methods
  • 5.4 Runge-Kutta Methods
  • 5.5 Multistep Methods
  • 5.7 Adams-Moulton Multistep Method
  • 5.9-5.10 Convergence and Errors
  • 6. BVPs
    • 6.1 Shooting Method
    • 6.2 Solution of Sets of Equations

1. Estimate the maximum step size h for Euler's Method applied to the first order ODE

   y'(x) = f(x, y(x))
   

   needed to make the absolute value of the theoretical global discretization error less than tol = 0.45e-2, given the (x,y)-dependent bounds

   |y| < 0.5223,  |f(x,y)| < 0.5649,
   
   |fx(x,y)| < 0.1432|x||y|,  
   |f_y(x,y)| < 0.04321|x|2,
   

   x on [-0.5710, +0.4710].

   ---

   Answer:
   
   Bounds: B=0.5649; B_x=0.04270; B_y=0.01408; B'=B_x+B*B_y=0.05065;
   so that h > tol*1*B_y/(B'*(exp(B_y*(b-a))-1)) > 0.1693 = h_max all using "=(r4)" precision.

   ---

   (S99a)

2. Using the 4th Order Runge-Kutta Method, numerically approximate the solution to the nonlinear IVP by

   y'(x) = 0.4710 y(x)(1 - 0.05444 y(x)/(1+x2)) ,  y(1) = 200 ,
   

   with h = 0.025 for two x-steps. Tabulate

   n		Xn	Yn		RK1	RK2	RK3	RK4		Xn+1	Yn+DeltaYn
   0
   1

   for all values calculated.

   ---

   Answer:

   n		Xn	Yn		RK1	RK2	RK3	RK4		Xn+1	Yn+DeltaYn
   0		1.000	200.0		-10.46	-9.714	-9.756	-9.703		1.025	190.2
   0		1.025	190.2		-9.069	-8.455	-8.487	-7.927		1.050	181.7

   ---

   (S99a)

3. Using the Algebraic BVP Method, numerically approximate the solution to the linear BVP,

   y''(x) + x y'(x) + 3 y(x) = -0.2 x2 , 
   
   y(1.25) = 2.5 & y(2.00) = 1.5 ,
   

   using finite central differencing of derivatives with h = 0.25. Tabulate

   n		0	1	2	3
   Xn
   Yn

   In solving, use the Thomas tridiagonal elimination algorithm. Sketch your approximation Y_n in the $xy$-plane.

   ---

   Answer:

   n		0	1	2	3
   Xn		1.250	1.500	1.750	2.000
   Yn		2.500	2.523	2.118	1.500

   The sketch of the approximation using linear interpolation is left as an further exercise.

   ---

   (S99a)

Caution: These old problems may have different answers due to different methods and chopping truncation used, so do not worry if your answers are a few least significant digits different.

For problems in the Gerald and Wheatley text,
• try Chapter 5 Exercises (pp. 455-466, 5th ed.):
  
  #5,6,7,8,10,11,12,13,14,20,29;
• try Chapter 6 Exercises (pp. 534-539, 5th ed.)
  
  #1 (use RK4, not RKF),3,8,9,13.
1. Consider the initial value problem:
   
    
   ODE: y'(t)=y(t)/(4+y(t)) = f(y(t))
    
   IC: y(0) = 1
    
   
        Let h =.5 and t(i) = i*h for i = 1 TO n. Find Y(i) = Y(t(i)) using 4-digit exam precision with the following methods
        (here, F(i) = f(Y(i)):
   
   (a) Euler's with n=3 tabulating

                  t(i), Y(i), F(i);
   

   
   (b) 4TH order Runge-kutta with n=2, tabulating

                   Y(i), RK1, RK2, RK3, RK4
   

   where needed.

   (Ans.:
   i  t(i)   YEuler(i)  FEuler(i)  YRK4(i) RK1(i) RK2(i) RK3(i) RK4(i)
   0  0.0    1.000      0.2000     1.000   0.1000 0.1039 0.1041 0.1081
   1  0.5    1.100      0.2156     1.104   0.1081 0.1122 0.1124 0.1165
   2  1.0    1.207      0.2318     1.216   ______ ______ ______ ______
   3  1.5    1.322      0.2484                   .......
   

2. Estimate the step size, h, needed to make the maximum global discretization error < 0.5e-3 of Euler's Method for t in (0,10) when
   
    
   |(d/dt)f(y(t),t)| < exp(t) & |partial(f)/partial(y)| < 0.5.
    

   (Ans.: ???)
   

3. Approximate the solution to the ODE:
   
    
   y'(t) = y(t)/(2+y(t)), y(0) = 1
    
   using Euler's Method for i = 0 to 6 stepping with h = 0.25. Tabulate your answers chopped to 4 significant digits with columns for

               i, t(i) =i*h, y(t(i)), y'(t(i)).
   

   (Ans.:   i  t(i)  y(t(i))  y'(t(i))
            0  0.00   1        0.3333
            1  0.25   1.083    0.3512
            2  0.50   1.170    0.3690
            3  0.75   1.262    0.3868
            4  1.00   1.358    0.4044
            5  1.25   1.459    0.4217
            6  1.50   1.564    0.4388  )
   

4. Find a numerical approximation to the solution of the initial value problem:
   
    
   y'(t) = f(t,y(t)) = exp(-t)/(1+y(t)), y(0) = 1,
    
   using a 4th order Runge-Kutta Method with h = 0.1 for t = 0 to 2*h. Tabulate

   	 t, Y, RK1, RK2, RK3, RK4
   

   at each step except the last, reporting only items needed. Use a modification of 4 digit precision.

   (Ans.:        t     y      RK1      RK2      RK3      RK4
                0.0  1.000  0.05000  0.04697  0.04700  0.04420
                0.1  1.047  0.04420  0.04159  0.04162  0.03919
                0.2  1.088                                                       )
   

5. a) Suppose that
   
    
   y'(t) = f(y(t),t) = exp(-y)/sqrt(1+t**2), y(0)) = 1
    
   and that the local chopping is bounded by 5.e-3, find the optimal step size h = 1/n for euler's method when 0 < t < 1. (Hint: Show that -1.5<(d/dt)f(y(t),t) < 0, y > 0, T > 0, getting a similar bound for f_y(y(t),t), and use this to estimate the Euler discretization error).

   
   b) For the ODE plus the IC in part a), illustrate Euler's Method with h = 0.1, tabulating

                 t,  y,  f(y,t) for t = 0 to 5*h.
   

   (Ans.:   a)  optimal h =0.8165e-2 (Caution:  Exponential error term not used)
            b)     t      y    f(y,t)
                  0.0   1.000  0.3678
                  0.1   1.036  0.3531
                  0.2   1.071  0.3360
                  0.3   1.104  0.3175
                  0.4   1.135  0.2984
                  0.5   1.164         )
   

6. Formulate an algorithm for a 4th order Runge-Kutta Method to numerically solve the initial value problem for autonomous systems of ordinary differential equations,
   
    
   y'(t) = f(y(t)), y(0) = C
    
   where y(t) = [y(i,t)]_{NX1} and f(y(t)) = [f(i,y(t)]_{NX1} are vectors with N components. Illustrate your algorithm by solving the system
   
    
   Y'(1,t) = Y(1,t)*(12-3*Y(1,t)-Y(2,t)), Y(1,0) = 1,
    
   Y'(2,t) = -Y(2,t)*(2-Y(1,t)), Y(2,0) = 1,
    
   with h = 0.1 for t = 0 to 1*h. Use a modified 4 digit precision.

   (Ans.:   ???)
   

7. Solve the following IVP numerically:
   
    
   y'(t) = sqrt(t)*y(t)**2, y(0) = 1,
    
   with h=0.1 and t(n)=n*h for n=1 to m using the following methods:
   
   a) First order Euler's Method with m = 4 tablulating

   
       t, Y,  F;
   

   
   b) 4th Order Runge-Kutta (RK4) Method for M=2, tabulating

   
           n, t, Y, RK1, RK2, RK3, RK4;
   

   (Ans.:
   n   t   YEuler   FEuler   YRK4   RK1     RK2      RK3      RK4
   0  0.0   1.000   0.0000   1.000  0.0000  0.02236  0.02286  0.03308
   1  0.1   1.000   0.3162   1.020  0.03290 0.04160  0.04195  0.05043
   2  0.2   1.031   0.4753   1.061    --      --       --       --
   3  0.3   1.078   0.6354
   4  0.4   1.141   0.8253                                            )
   

8. Estimate the step size h needed to make the absolute value of the global theoretical discretization error on (0,3) less than 1.e-4 when
   
    
   -3*t < f_x(t,y(t)) < 2*t
    
   -3*t < f_y(t,y(t)) < 2*t
    
   abs(f(t,y(t))) < 1.
    

   (Ans.:   abs(f') ≤ abs(f_x)+abs(f)*abs(f_y) < 9+1*9=18;   Using
            abs(global error) = E_n < m*h*(exp(K*(x_n-x_0))-1)/(2*K),
   	 then h <(4C) 1.879E-16)
   

9. Consider the BVP
   
    
   y" + y' + x*y =0 , y(0) = 0, y(5) = 1.
    
   Apply the central difference BVP Method with h=1 to convert the BVP into a linear algebra problem with the proper numerical coefficients, but do not solve the matrix problem.

   (Ans.:
       ( -1.0   1.5   0.0   0.0 ) ( Y(1) )      (  0.0 )
       (  0.5   0.0   1.5   0.0 )*( Y(2) )   =  (  0.0 )
       (  0.0   0.5   1.0   1.5 ) ( Y(3) )      (  0.0 )
       (  0.0   0.0   0.5   2.0 ) ( Y(4) )      ( -1.5 )     )
   

10. Find a numerical approximation to the solution of
    
     
    y'(t) = f(t,y(t)) = y(t)*exp(y(t)), y(0) = -1,
     
    using a 4th order Adams' (-Bashforth) Predictor-Corrector Method with a single derivative function evaluation following each "P" or "C" step, h = 0.25 for t=i*h and i = 1 to 5. Use just enough Euler steps to start the Adams' Method. Tabulate both y and f at both predictor and corrector steps.

    (Ans.:   i  method  t(i)  y(i)   f(i)
             0   I. C.   0   -1.000  -0.3678
             1  EULER   0.25 -1.091  -0.3664
             2    "     0.50 -1.182  -0.3624
             3    "     0.75 -1.272  -0.3565
             4    P     1.00 -1.360  -0.3490
             "    C       "  -1.360  -0.3490
             5    P     1.25 -1.446  -0.3405
             "    C       "  -1.446  -0.3405 )
    

11. (Master's Exam Spring 1998, NA#11) Estimate the step size  h  for Euler's Method needed to make the absolute value of the theoretical global discretization error less that 0.25e-1, given the bounds

    |f(x,y)| < 1.387,
    |f_x(x,y)| < 0.6795*|x|,
    |f_y(x,y)| < 0.5432*|x|,
    

    on [-2.25,1.35].

    {Optional Hint: It is permissible to exponentially approximate  (1+h*c)^n  for some bounded constant  c , where  n  is the number of steps.}

    Use 4 digit exam precision: chop recorded intermediate results only to 4 significant decimal digits and continue calculations with these results.

    ---

    Answer: The theoretical global error is such that

    |E| < |(B'*h/(2*By))*[exp((x_n-x_0)*By)-1]| < tol,
    

    where tol=0.25e-1; and the bounds are given by

    By=0.5432*2.25=(4ch) 1.222; 
    
    
    Bx=0.6795*2.25=(4ch) 1.528; 
    
    
    B=1.387; 
    

    and where B'=Bx+B*By=3.222. Solving the error inequality for the best step size h yields

    h < tol*2*By/[B'*[exp(By*(1.35-(-2.25)))-1]]=(4ch) 0.2359e-3 (answer).
    

    {Note: The less than "<" symbol used in the HTML form of this problem should really be the less than or equal "<=" symbol.}

    ---