Problem of the Day

**Aug 26** Assume
f and its first five derivatives are continuous. Using Taylor series one can
show that

f(a+h)-f(a-h) = 2hf'(a) + (h^3/3)f'''(a) + O(h^5)

where the "O" terms means stuff left over of the form (h^5)*(modest multiple of f's fifth derivative evaluated somewhere.) Rewrite the formula with h replaced by 2h. Now manipulate the two equations to get a recipe of the form

f'(a) = ( linear combo of f(a-2h), f(a-h), f(a+h), f(a+ 2h)) + O(h^4)

Implement a function of the form yp = Derivative4(f,a,h) that uses your recipe to approximate f'(a). Solution

--------------------

**Aug 28** What is the largest power of
10 that can be represented exactly in IEEE floating point? Harder: What is the
largest value of n so that n! can be represented exactly in IEEE floating point?
Write scripts that answer these questions. Hint

--------------------

**Sep 2** (a) Write a function
PlotDerPoly(x,y) that plots the derivative of the polynomial interpolant of the
data (x(i),y(i)), i=1:n. Assume

that x(1)<...< x(n) and that the plot is across [x(1),x(2)]. Make effective use of the Matlab built-in functions polyfit and polyval. (b) Write a function a = InterpV(x,y,u,v) where x and y are column n-vectors, u and v are scalers, and a is a column n-vector so that if

p(z) = a(1) + a(2)(z-u)/v + a(3)*((z-u)/v)^2 + ... + a(n)*((z-u)/v)^(n-1)

then p(x(i)) = y(i), i=1:n. Assume that the x(i) are distinct. Hint

--------------------

**Sep 4** Play with
ShowCheby. Then Modify
ShowRungePhenom so that it
displays the interpolating polynomial based on the Chebychev points. In AG, go
to pages 323-324 and know the answers to 0(a,b,c,d,e,i,j,k,m,o,r,t)

--------------------

**Sep 9 **Assume that x is a column n-vector
with x(1) < ... < x(n). (a) Write a fully vectorized function w = BaryWeights(x)
that returns a column n-vector w with the property that w(j) = (x(j) - x(1))
*...* (x(j) - x(j-1))*(x(j)-x(j+1)) *...*(x(j) - x(n)) for j=1:n. (b) Write a
fully vectorized

function zVals = LagrangeEval(u,w,x,y)

that accepts column n-vectors w and x and column m-vector u and returns a column m-vector zValswith the property that zVals(i) = p(u(i)) where

p(z) = ( y(1)*w(1)/(z-x(1)) + ... + y(n)*w(n)/(z-x(n) ) / (w(1)/(z-x(1)) + ... + w(n)/(z-x(n) ) )

Assume that x(1)< ... < x(n) and no u-component is an x-component. See page 305 in Ascher and Greif. Hint

--------------------

**Sep 11** Assume that x and y are
column n-vectors and x(1)<x(2)<...<x(n). Suppose in Matlab we execute S =
spline(x,y). Let s(x) be the spline that is produced. Show how to compute the
integral of [s"(x)]^2 from x(1) to x(n). Hint

--------------------

**Sep 16 **Suppose S = spline(x,y). How would
you compute the biggest jump in the 3rd derivative of the spline
represented in S? Hint

--------------------

**Sep 18** Consider
pwLadapt. Suppose it is called
with hmin =0. And suppose you know that the |f''(x)| <= M2 for all x. pwLadapt
will return an x-vector consisting of the breakpoints. Give a lower bound for
min |x(k+1) - x(k)|. Hint

--------------------

**Sep 23** (a) Refer to Bstar. Let tau be the
integral of B_{*}(z) from -2 to 2. What is tau? (b) Refer to equation (2) n the
A3 handout. Given tau, what is the integral of s(z) from x(1) to x(n)? Assume
that the alpha(k) are known. Hint

--------------------

**Sep 25** Suppose [a,b] is given and
that we also have a vector alfa(1:n+1). Define the egree n polynomial p(x)
in terms of the Chebychev polynomials T_{0},T_{1},...,T_{n}as follows:

p(x) = alfa(1)*T_{0}(-1 + 2(x-a)/(b-a)) + alfa(2)*T_{1}(-1 + 2(x-a)/(b-a)) + ... + alfa(n+1)*T_{n}(-1 + 2(x-a)/(b-a))

How would you compute p'(x0) where x0 is given? Hint

--------------------

**Sep 30 ** The nxn dft matrix F is
defined by [F]_{kj} = omega^(k*j) where omega = cos(2pi/n) - i*sin(2pi/n).
(Subscripting from zero.) Show that if Fbar = conj(F), then Fbar*F = n*I.
Hint

--------------------

**Oct 2** Rewrite
CompQNC so that it involves a single
function call. Hint

--------------------

**Oct 7** Download
CompSimpson,
AdaptSimpson, and
ShowSimpson. (a) In AdaptSimpson,
why do we have Q = Q2 +errInQ2 and not just Q = Q2? (b) Modify
AdaptSimpson so that useful f-evaluations that are ``on tap'' can be passed in
the recursive calls. Hint

--------------------

**Oct 9** Write a function Q = CompGL(f,a,b,n)
that estimates the integral of f from a to b by applying the 2-point
Gauss-Legendre rule to each subinterval defined by z = linspace(a,b,n+1). Give
an error bound. Hint

--------------------

**Oct 16** The Matlab quadrature routine
quadl is exact for polynomials of degree m or less. Write a script that
determines m experimentally. Hint

---------------------------------------------------------------Mid term cut-off -------------------------------------------------------

**Oct 28 **In
ShowRabbits, replace the recipe for dr/dt with dr/dt = 2(1 - r(t)/R)*r(t) -
alpha*r(t)*f(t). Play with the constant parameter R and summarize what you
observe. Hint

--------------------

**Oct 30** How would you use ode45 so solve
(d^2 x /dt^2) + A*(dx/dt) + Bx = 0, x(0) = x0, (dx(0)/dt) = xp0 where A and B
are nxn matrices and x0 and xp0 are column n-vectors.
Hint

--------------------

**Nov 4** Download
ShowEvent. (a) Write and test an event
function that is triggered when the rabbit and fox densities are the same.
Hint (b)
Download ShowMarkov. Let E(t) = expm(Q*t).
Augment ShowMarkov so that it plots the first column of E(t) across the interval
[0,4] Hint

--------------------

**Nov 6** Consider dy/dt = A*y(t) with y(0) =
y0 and the AB2 method with stepsize h. If y_{k} and y_{k-1}
are exact, what can you say about y_{k+1}? Hint

--------------------

**Nov 11** For a fixed point iteration

z_{k+1} = G(z_{k}

to converge, we need |G'(z)|<=c < 1. Now consider using the fixed point iteration method for computing y_{k+1} via the AM1, AM2, AM3, AM4 and AM5 methods. What is the "G-function" in each case and give a condition (involving h) that guarantees convergence, Hint

--------------------

**Nov 13** Refer to the lecture code
ShowFlame. "This problem is only stiff in part of the interval [0,2/delta]".
Explain. What is the steady state solution if

y'(t) = alpha*y^2 - beta^2 * y^3 y(0) = delta

where alpha>0 and beta > 0. In th eoriginal problem, big things happen near t=1/delta. What about for the modified problem? Hint

--------------------

**Nov 18 ** Take a look at ShowTwoPtBVP.
Modify the subfunction TwoPtBVP(a,ua,b,ub,n,p,q,r) so that ua and ub prescribe
endslopes, i.e., ua = u'(a) and ub = u'(b). The additional linear equations
(u(2)-u(1))/h = ua and (u(n)-u(n-1))/h = ub mean that the linear system is now
n-by-n. Hint

--------------------

**Nov 20 ** Let T be the nxn tridiagonal
matrix with 2's down the diagonal and -1's down the subdiagonal and
superdiagonal. This matrix comes up in boundary value problems a lot. Suppose Tu
= f where is an n-vector of the form [ua ; zeros(n-2,1) ; ub] where ua and ub
are positive. Show that either u(1) or u(n) is the largest entry in the solution
vector. Hint

--------------------

**Nov 25 ** Consider the Rabbit-Fox
simulation where the initial fox distrinution is fox(0) = 10. How would you
determine the initial rabbit density rab(0) so that at a specifed time T, rab(T)
= 200? Just outlinet the overall method. Hint

--------------------

**Dec 2** Download FastPoisson and modify it
so that it solves u_xx + u_yy = 0 on the rectangle but that u(x,y) agrees with
f(x,y) on the boundary. Hint

--------------------

**Dec 4** On page 531 of Ascher&Greif,
problem 0, (A)-(L),(N)

--------------------