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)
--------------------