MCS 471 Lecture Forty Four

(Partial) Answers to Review Questions of Chapters 6 and 7

You may use a formula sheet: version in pdf format .

Apply the method of underdetermined coefficients to derive an Adams-Moulton formula which uses three function evaluations. What is the order of the local error (in terms of the step size h) of the method you just derived? Justify your answer.

Answers:

If we solve dy/dx = f, then we can integrate dy = f dx over the interval [x(n),x(n)+h]. At the left of the equality we obtain y(n+1) - y(n), at the right we set up a quadrature formula at the points x(n)-h, x(n) and x(n) + h. To simplify the calculations, we may set x(n) = 0. Requiring an algebraic degree of precision of two, we find the weights w1, w2, w3 in the quadrature formula as the solution to this linear system:

             h
             /
            |   1 dx =  h  = w1 + w2 + w3
            /
             0

             h           2
             /          h  
            |   x dx = --- = w1 (-h) + w2 0 + w3 h
            /           2
             0

             h           3
             /   2      h           2             2
            |   x dx = --- = w1 (-h) + w2 0 + w3 h
            /           3
             0

The solution

             -1              +8              +5
      w1 =  ---- h     w2 = ---- h     w3 = ---- h
             12              12              12

defines the three point Adams-Moulton method:

              h  (     (            )       (        )       (            ) )
  y   = y  + --- ( - f ( x   , y    ) + 8 f ( x , y  ) + 5 f ( x   , y    ) )
   n+1   n    12 (     (  n-1   n-1 )       (  n   n )       (  n+1   n+1 ) )

The interpolating polynomial p(x) has the following error (for x_n = 0):

        (4)
       f   (a)
      --------- (x - h) ( x - 0 ) ( x - (-h) ),  for some a in [-h,+h].
         4!

Integrating this error over [0,+h] gives

        (4)      h                 (4)         4        2
       f   (a)  /   3    2         f   (a) (  h     2  h  )      4
      --------- |  x  - h x dx = --------- ( --- - h  --- ) = O(h )
         4!     /                   4!     (  4        2  )
                 0

Thus we derived a method whose local error is O(h^4) and we can write

                      h  (                        )      4 
    y(x   ) = y  +  ---- ( - f    + 8 f  + 5 f    ) + O(h )
       n+1     n     12  (    n-1      n      n+1 )

The global error is O(h^3) and therefore we speak of a third-order method.

Consider the equation y" = x+y, with y(0) = 1 and y'(0) = 2. Apply a predictor-corrector method, using three points in each step to approximate y(0.5), with step size h = 0.1.
Hint:
Note that we must apply a third-order Runge-Kutta method to get started. This problem is analoguous to what we did in Lecture 37, except for that f is now a system: for every x you have y(x) and y'(x).
Derive a bound on the step size h so that Euler's method converges on the problem y' = lambda y, y(0) = 1, for lambda < 0. Can you also derive a bound on h for the modifield Euler's method?
Answers:
Euler's method on this problem amounts to y(n+1) = (1+lambda*h)*y(n). In order for the method to converge for negative lambda, y(n) must go to zero, or equivalently |1 + lambda*h| < 1. This leads to the upper bound 2/(-lambda) for h.
For the modified Euler's method, we find
```
               (1+lambda*h/2)
     y(n+1) = ---------------- y(n)
               (1-lambda*h/2)
```
Here we derive no constraint on h, other than the trivial condition h>0.
Consider the following boundary-value problem:
```
      y" - 2 x y = x^3,     y(0) = 1, y(1) = -1.
```
With the shooting method we translate this problem into initial-value problems.
1. With our first guess y'(0) = 1 we find 3.12767 at x=1. Our second guess y'(0) = -1 yields -0.320365 at x=1. What is your next guess for y'(0) in the shooting method?
2. Suppose we modify the differential equation above into y" - 2 y^2 = x^3. Is the modified problem easier or harder than the original one? Justify your answer.
3. Suppose, instead of the conditions y(0) = 1 and y(1) = -1 we have mixed boundary conditions: y'(0) = 1 and y(1) = -1 for the original ODE above. Describe, in general (but in detail), how the shooting method works to solve this problem.
Answers:
1. We solve a linear interpolation problem:
```
   p(1) = 3.12767    p(-1) = -0.320365

           z - 1                 z + 1
   p(z) = ------- (-0.320365) + ------- 3.12767
          -1 - 1                 1 + 1

        = 1.72402 z + 1.40365

   p(z) = -1   => z = -1.39422
```
2. The modified differential equation is nonlinear. Therefore, the linear interpolation will no longer yield the correct value for y'(0). Thus the modified problem is harder than the original one.
3. If we are not given values for y(0), then we make two guesses, we call A and B our two guesses for y(0). With y(0) = A, y'(0) = 1 we solve an initial value problem and obtain as solution yA. Similarly, we obtain yB as solution to the initial value problem with initial conditions y(0) = B, y'(0) = 1. Let p be the interpolating polynomial through (A,yA(1)) and (B,yB(1)). Solving p(x) = -1 gives a new guess for y(0).

Consider the boundary-value problem

      y" + x y' = x^2,     y(0) = 1, y(1) = -1.

With finite-difference approximations we translate this problem into a linear system.

Give the linear system one has to solve for h = 0.2.
Suppose, instead of the given boundary conditions, we are given y'(0) = 1 and y'(1) = 1. Given the linear system one has to solve for h=0.2.

Partial Answers:

We use the following approximations for the derivatives at the i-th point:

              1   (                    )            1   (            )
       y" = ----- ( y    - 2 y  + y    )     y' = ----- ( y   - y    )
        i    h^2  (  i+1      i    i-1 )      i    2 h  (  i+1   i-1 )

The the equation above becomes (after multiplication with h^2):

                             h  (             )    2  2
    y    - 2 y  + y    + x  --- ( y    - y    ) = x  h 
     i+1      i    i-1    i  2  (  i+1    i-1 )    i

or equivalently:

   (         h  )               (         h  )         2  2
   ( 1 + x  --- ) y    - 2 y  + ( 1 - x  --- ) y    = x  h
   (      i  2  )  i+1      i   (      i  2  )  i-1    i

With the values for y given at the boundary, we have four unknown values for y, we call them y1, y2, y3, y4, with corresponding values for x at x1 = 0.2, x2 = 0.4, x3 = x.6, and x4 = 0.8. The linear system we have to solve is the following:

 [    -2     1+0.2*0.1     0        0     ][ y1 ]        [ 0.2^2 - (1-0.2*0.1)*1  ]
 [                                        ][    ]        [                        ]
 [ 1-0.4*0.1    -2     1+0.4*01     0     ][ y2 ]        [ 0.4^2                  ]
 [                                        ][    ] = .2^2 [                        ]
 [     0     1-0.6*0.1    -2    1+0.6*0.1 ][ y3 ]        [ 0.6^2                  ]
 [                                        ][    ]        [                        ]
 [     0         0     1-0.8*01    -2     ][ y4 ]        [ 0.8^2 - (1+0.8*0.1)(-1)]

For derivative boundary conditions, we introduce fictitious nodes at x(-1) and x6. We eliminate the values at those nodes, applying central difference formulas for the derivative and using the derivative conditions on the boundary. The system we have to solve is now of dimension six instead of four, with unknowns y0,y1,y2,y3,y4,y5. If we set up the equations (like in the first part for this question), we obtain for i=0, the point y(-1) and, for i=5, the point y6. To eliminate those conditions, we apply the boundary conditions, using central finite-difference approximations for the derivative:
```
            y(0+h) - y(0-h)     y1 - y(-1)
   y'(0) = ----------------- = ------------ 
                  2h               2h
```
so we find y(-1) = y1 - 2*h*y'(0).
```
            y(1+h) - y(1-h)     y6 - y4
   y'(1) = ----------------- = --------- 
                  2h               2h
```
so we find y6 = y4 + 2*h*y'(1).

The equation below depends on a parameter k:

     y" - 3 y' + 2k^2 y = 0,     y(0) = 0, y(1) = 0.

With finite-difference approximations we translate this problem into an eigenvalue problem.

Give the eigenvalue problem one has to solve for h = 0.2.
Explain how to approximate the eigenvector whose corresponding eigenvalue lies closest to 1.

Answers:

We use the following approximations for the derivatives at the i-th point:

              1   (                    )            1   (            )
       y" = ----- ( y    - 2 y  + y    )     y' = ----- ( y   - y    )
        i    h^2  (  i+1      i    i-1 )      i    2 h  (  i+1   i-1 )

The the equation above becomes (after multiplication with h^2):

       y    - 2 y  + y    - 3/2 h y    + 3/2 h y    + 2 k^2 h^2 y  = 0
        i+1      i    i-1          i+1          i-1              i

or equivalently:

  (3h/2-1) y    + 2 y  - (3h/2+1) y    =  2 k^2 h^2 y 
            i+1      i             i-1               i

With h = 0.2, we have four internal points inside [0,1], we call them y1,y2,y3,and y3. Thus we must solve

    [    2    3h/2-1     0       0   ] [ y1 ]             [ y1 ]
    [                                ] [    ]             [    ]
    [ -3h/2-1     2   3h/2-1     0   ] [ y2 ]             [ y2 ]
    [                                ] [    ] = 2 k^2 h^2 [    ]
    [   0    -3h/2-1     2    3h/2-1 ] [ y3 ]             [ y3 ]
    [                                ] [    ]             [    ]
    [   0       0    -3h/2-1     2   ] [ y4 ]             [ y4 ]

Setting lambda = 2 k^2 h^2, the system above is an eigenvalue problem which defines four eigenvalues and eigenvectors.

To find the eigenvector whose corresponding eigenvalue lies closest to 1, we apply the shifted power method. If the matrix is denoted by A, then we execute the following iteration, starting at some vector x(0), and continueing with
```
           (A - 1 I) x(k+1) = x(k)
           x(k+1) = x(k+1)/norm(x(k+1))
```
where I is the identity matrix.

FINAL EXAM is in BH 0309 on Monday 29 April 2002, from 1 till 3PM.