-------------------------------------------- Math. 505 Fall'12 - Preliminary Syllabus -------------------------------------------- updated: 11/05/12 problems with a single number are the same in eds. 1,2 problems listed as {m/n} are m in ed. 1 and n in ed. 2 -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- Week 1: Introduction: machine arithmetic, errors and conditioning. Set 1: (Ch.1) 2,3,5,6,{40/42};programming: 2,6,{14/17} Due: 9/6/12 {ed1/ed2} 1. (Matlab)>> help eps EPS Spacing of floating point numbers. D = EPS(X), is the positive distance from ABS(X) to the next larger in magnitude floating point number of the same precision as X. X may be either double precision or single precision. For all X, EPS(X) is equal to EPS(ABS(X)). EPS, with no arguments, is the distance from 1.0 to the next larger double precision number, that is EPS with no arguments returns 2^(-52).... HINT: for machine problem {14/17} it is possible to get a closed-form solution for this 2-term recurrence relation. First look for a "homogeneous" solution, i.e. a solution of the equation Y(n) = n*Y(n-1) then find a particular solution by method of undetermined coeffs: let Z(n) = C(n)*Y(n) then the general solution is y(n) = Z(n) + c*Y(n) where c is an arbitrary constant. Important thing to remember is that even if you think you are computing Z(n) which is bounded, there is numerical contamination, adding a small part of Y(n) which blows -------------------------------------------------------------------------------- Week 2-6: Approximation of functions: L2 vs. Interpolation. Chebyshev and Lagrange interpolation; splines. Convergence and error minimization. Set 2: (Ch.2) {7/8}, {13/14}, {24/27}, {56/62}, {58/64} ({1st Ed/2nd Ed.}) machine 7,{9/8}. Due 10/9/12 -------------------------------------------------------------------------------- Week 6-8: Numerical differentiation and integration; the connection to the interpolation problem and extrapolation. Set 3: (Ch.3) {6/7},{7/8},{12/13},{22/23},{35/36},{36/37} machine: posted part A (see mach.prob.6 for ideas); {3/4} NOTE: the answers given by Gautschi, 1st Ed., for machine problem 4 are only single precision accurate. The double precision answers are: (corrected in 2nd Ed.). \Int_0^1 cos{x}/sqrt{x} dx = 1.80904847580054 \Int_0^1 sin{x}/sqrt{x} dx = 0.62053660344676; Due: 10/25/12 Sample output for Rhomberg integration: f(x) = x^2*exp(-x); 0<=x<=1 h = 2^(-k), k = 0,1,2,3,...9 k Trapezoidal Rhomberg Exact 0 0.18393972 0.18393972 0.16060279 1 0.16778619 0.16240168 0.16060279 2 0.16248841 0.16061053 0.16060279 3 0.16107990 0.16060280 0.16060279 4 0.16072243 0.16060279 0.16060279 5 0.16063272 0.16060279 0.16060279 6 0.16061028 0.16060279 0.16060279 7 0.16060467 0.16060279 0.16060279 8 0.16060326 0.16060279 0.16060279 9 0.16060291 0.16060279 0.16060279 -------------------------------------------------------------------------------- Week 9-10: Nonlinear Equations: bisection & Sturm sequences, secant and Newton's methods, fixed point iterations. Set 4: 11, 26, 35, 36; machine 2, 4, 5. Due Thu. Nov. 15 Week 11-15: The solution of ordinary differential equations (ODEs) 11-13: Initial value problems: implicit vs. explicit, single and Set 5: 7, 12, 18, 19; machine 3. Due Thu. Nov. 29 performance of RK3, RK4 and RK5 on the 4 cases (shown are the number of stepsize halvings to achieve tolerance, h and maximum absolute error (RK4 and RK5 results are shown for comparison) RK3 iter h max_err 06 0.00312500000000 0.00000009625056 09 0.00039062500000 0.00000006461166 09 0.00039062500000 0.00000032468692 10 0.00019531250000 0.00000011546136 RK4 iter h max_err 05 0.00625000000000 0.00000004225998 07 0.00156250000000 0.00000032052294 08 0.00078125000000 0.00000011366661 08 0.00078125000000 0.00000023098397 RK5 iter h max_err 03 0.02500000000000 0.00000006660244 06 0.00312500000000 0.00000001659540 06 0.00312500000000 0.00000015089578 07 0.00156250000000 0.00000001934639 13-14 multistep methods; stability and stiffness. Set 6: (see posted) due by noon, Monday, Dec. 10 You must submit an electronic copy of your matlab code so we can reproduce your calculations. ----------------- EXTENSION: Monday, December 10, 19:30 *** Electronic Submission Required!!! *** I will accept late assignments (or corrections/additions to what you have already turned in) submitted electronically before midnight, Friday, Dec. 14. ----------------- Also, note that due to the perturbed moon/earth mass ratio, if you decide to hunt for alternative periodic solutions you will need to employ an iterative scheme, to search for alternative values for the period T and the initial condition for y'. This part is not required! (Hint: be sure you verify that the main part of your scheme works correctly before you focus on the adaptive part. Depending on the scheme you choose - for my own calculation I used a 6th order Milne scheme: Y7 = y1 + h*(a*f2+b*f3+c*f4+b*f5+a*f6) G7 = F(x7,Y7) y7 = y3 + h*(A*f3+B*f4+C*f5+B*f6+A*G7) a = 33/10; b = -21/5 ; c = 39/5 ; A = 14/45; B = 64/45; C = 24/45; combined with a 5th order RK scheme for initialization and step-halving. What happens to your calculations after a period without adaptivity? Mine goes unstable unless I use a VERY small time step!) As for the error tolerance (||U-U*||) it is hard to estimate the value here since the potential function has a singularity (although there are ways around it-using, say conservation of energy to estimate the closest the trajectory may come to one of the two attracting centers. Much easier to use an absolute bound, experimenting with the value until you compute an orbit that closes and remains closed after 4 revolutions. I suggest starting with a modest value, like 10E-3 and decrease if necessary. a step control strategy: assuming you are using scheme AB of order p; need p steps to get p+1 step. Must store 2*P-2 values to accomodate step doubling. define tolerance eps, stepsize h. startup: given Y1, compute F1 -take p-1 steps with RK scheme. After each step k, let Y(k+1) = U, F(k+1)=YP(U). -take 1 step with AB. Find U, U*. Compute S = ||U-U*|| -If S > eps, let h = h/2 and repeat from begining, until a stepsize is found that allows enough accuracy. -If S < .1*eps then the stepsize can be increased. However for that you need some additional past values: a P+1 (P+1-a)*H +----+----+----+----+----+----+----+----+----+----+ -3 -2 -1 0 1 2 3 4 5 6 7 +---------+---------+---------+---------+---------+ 1 2 3 4 5 6 1 P (P-1)*2*H gives a = 3-P =-3 for P=6 i.e. for a scheme of order 6 you need to keep the 6 values/derivatives needed to advance to the next step, but also 4 additional earlier values in case you need to double. Obviously after a doubling (as after a halving) you must take at least P-2 additional steps to store enough past info before doubling again. After a step is taken, if result is acceptable, cycle values back: Y(k+1) -> Y(k); U -> Y(P). I find it easier to index backward (i.e. Y(1) is latest value, Y(2) previous, etc.; for my 6th order scheme, I needed two 10x4 arrays, for previous values and derivatives. 14-15: Boundary value problems: shooting or Ax=b (or F(x)=b)? -------------------------------------------------------------------------------- --------------------------------------------------------------------------------