Final Exam Date:
Numeric (Matlab) portion: May 2 (Mon) 1:00-2:30pm, PC Lab 36 (usual meeting time/place)
Analytic portion: May 5 (Thu) 3:15-5:45pm, Kupf 108


Final Exam - Theory Portion (May 5, 2005): [ exam ]    Solution: [ 1 ]   [ 2 ]   [ 3 ]

Final Exam - Matlab Part
1 stephensen.m Matlab function stephensen(f, x0, tol) calculates the root of equation f(x)=0, using initial condition x0. execute stephensen( inline('exp(-x)-x'), 1, 1e-8 );
2 dice.m Matlab function dice(diceSum, nTrials) calculates the average number of two-dice rolls required to get a sum of diceSum, using a Monte-Carlo simulation with nTrials trials. Execute dice(4, 100), dice(4, 1000) and dice(4, 1000). Exact result is 12.
3a,b rule3_8.m Function rule3_8(f, a, b, N) calculates the integral of inline function f from a to b using the 3/8 rule with N subintervals.
problem3.m Matlab code for problem 3. The exact slope of the line is -1.5, indicating that the error decreases as N-3/2. Note that the convergence to exact solution would be much faster if the integrand did not have an infinite derivative at x=0. In general, the error of the 3/8 rule decreases as N-5 (faster than for the Simpson's rule)



Midterm (March 9, 2005): [ exam ]    Solution: [ 1 ]   [ 2 ]   [ 3 ]   [ factorial m-function ]  


Homework Solutions: [ 3 ] [ 4 ] [ 5 ] [ 6 ] [ 7 ] [ 8 ] [ 9 ] [ 10 ] [ 11 ] [ 12 ]


Go to Lecture Notes


Homework #3
1 hw3_1.m Limit of (ex-1) / x for x -> 0
2 hw3_2.m Limit of (root(x2+1) - 1) / x2 for x -> 0
3 IntCosX2.m Matlab function integrating cos(x2) with a tolerance of 10-5
4 IntLn.m Matlab function integrating ln(1+x) with a tolerance of 10-5
2, 3d HW3_2_3d.doc Analytic results for problems 2 and 3d
Homework #4
1a hw4_1.m Table comparing f(x)=.5 x2 ((x-1)-1-(x+1)-1) and g(x)=x2 (x2 - 1)-1
1b g(x) is clearly more accurate - it gives a correct limit of 1 as x -> Inf, and remains greater than 1 (see table). f(x) gives poor results since it contains a difference of two decreasing values. Even a small round-off error becomes significant as this difference gets smaller and smaller.
2a binomial.m Matlab function calculating explicitly the binomial expansion (x + y)n
2b hw4_2b.m Table of (x+y)15 values calculated using binomial.m, for x=1:20, and y=1-x
For n=15, the binomial expansion only works when x < 7, since it involves powers of x and y up to 15, which are very large for x ≥ 7. For instance when x=7, the k=7 term is larger than 1015 ( |a7| = x8 y7 15! / (7! 8!) = 1.04 * 1016 ), resulting in a loss of significance of order one, so the error is comparable to the total sum.
2c hw4_2c.m Table of (x+y)n values calculated using binomial.m, for n=1:20, x=10, and y=-9
For x=10 and y=-9, the binomial expansion only works when n < 14, as we would expect, since loss of significance becomes comparable to one when terms become larger than 1015. When n=14, the k=7 term is |a7| = 107 97 14! / (7! 7!) = 1.64 * 1017
3 mySin.m Matlab function calculating sin(x) by summing the Taylor series
hw4_3.m Table comparing mySin(x) to built-in sin(x) values, for x=0:0.1:pi. Note the almost perfect agreement between the two columns
hw4_3b.m Comparison for larger angles: mySin(x) fails for x around 35-40 radian. This is because the terms in the Taylor series are becoming very large as x increases (larger than 1015), leading to round-off errors comparable in value to the sum itself
Homework #5
1a hw5_1a.m Table of trigonometric error sin(1+πn)+sin(1), where n is an odd integer.
1b Equation sin(x + πn)+sin(x)=0 will fail miserably when all information about x in sin(x + πn) is lost. This happens when the round-off error in representing πn becomes comparable to x. The round-off error of πn is less than 10-15πn. Therefore, the problem occurs when x < 10-15πn, which for x=1 is satisfied when n > 1e15 (we don't have to be too accurate here, since we only know the upper bound on error, and not its exact magnitude). Let's check using Matlab: calculation sin(1+pi*(1e15+1))+sin(1) yields -0.06, which may seem tolerable, but for a slightly higher n=2e15+1 we get sin(1+pi*(2e15+1))+sin(1)=0.47, which definitely represents complete failure.
1c newSin.m Matlab function calculating the sine using its Taylor series. The argument is normalized to lie in the interval (-π, π) before the series is summed, improving the results for large values of the angle.
hw5_1c.m Table of trigonometric error sin(1+πn)+sin(1), where n=10m+1, and m=1:17, for the built-in sine and the newSin functions. This program also plots the logarithm of this error vs. m (which itself is the log10 of n). The complete failure occurs when log(error) is close to zero, since this corresponds to error of order one. On the graph this clearly corresponds to m of about 15-16.
2a hw5_2a.doc Proof of limit (1+x/N+(x/N)2/2)N -> exp(x) as x -> Infinity
2b limExp.m Matlab function calculating the exponential function using the limit calculation. Note the funny header: this function will return two numbers - the value of the exponential function, and the highest value of m used to find it.
2c hw5_2c.m Comparison of limExp(x) to the built-in function exp(x), for x from 1 to 20. Note the interestingly sudden failure at x=18.
The explanation for failure at x=18 is simple: as x increases, higher and higher N is needed to calculate the limit (see how m grows with increasing x). In fact, we know when the expression under the limit will fail completely: this will happen when x/N will be of order 1e-15 to 1e-16, since we know that such a small number will get completely rounded off when added to 1. Now, for x between 10 and 20, this occurs when N > 1e16 or 1e17, which corresponds to 2^53 to 2^56. From the table you can see that this is exactly the region of m values where the method failed!
Homework #6
1 hw6_1.m Solution is yn=(2/3)2n
2 hw6_2.m Solution yields r1=1/2, r2=3/2, c1=2, c2=0, so yn=2/2n=1/2n-1
3 hw6_3.m Constant solution of non-homogeneous equation is C = -1, so yn=c1(1/2)n+c2(3/2)n-1; initial conditions yield c1=5/2, c2=1/2, so yn = (5/2)(1/2)n + (1/2)(3/2)n - 1 = (5 + 3n) / 2n+1 - 1
2, 3 The difference in problems 2 and 3 is zero since there are only factors of 2 in the denominator - therefore, in this case there are no round-off errors (at least untill yn become larger than 253)
Homework #7
Homework #8
1a hw8_1a0.m yn+1 = 64 yn / (1 + yn)3. Equilibria: yeq=0 and yeq=3. Stability: f '(y) = 64 (1 - 2 y) / (1 + y)4
yeq=0: f '(0) = 64
(monotonically unstable), yeq=3: f '(3) = -5/4 (oscillatorily unstable);
hw8_1a3.m
1b hw8_1b.m yn+1 = yn ln(yn). Equilibria: yeq=0 and yeq = e = 2.718....
Stability of non-zero equilibrium: f '(y) = ln(y) + 1 -> f '(e) = 1 + 1 = 2 (monotonically unstable)
2 yn+1 = yn er (1 - yn). Equilibria: yeq=0 and yeq = 1. Stability: f '(y) = (1 - r y) er (1 - yn)
Zero equilibrium: f '(0) = er: stable if r < 0; Non-zero equilibrium: f '(1) = 1 - r: stable if 0 < r < 2
hw8_2.m Matlab portion of problem 2 (pick two values from just inside the interval 0 < r < 2, and two values of r from just outside this interval)
Homework #9
1, 3 plotRicker.m Matlab function plotRicker(r, y0, n, s) directly computes and plots the solution of the Ricker's equation with parameter r, initial condition y0=y0, for n steps, using line type s.
1 hw9_1.m Plots the directly computed solution to the Ricker's equation for r=2.2, r=2.4, and r=3.4.
2 hw9_2.m This Matlab function calculates the bifurcation diagram for the Ricker equation yn+1=yn er (1 - yn)
hw9_2.jpg The bifurcation diagram for the Ricker equation; note the similarity with the logistic equation bifurcation diagram
3 hw9_3.m The "butterfly effect": note how the slight difference in initial condition is amplified as n grows.
4 butterfly.m Matlab function butterfly(y0, d, r) returns the iteration number at which the solutions of the Ricker's equation with parameter r and initial conditions y0 and y0+d deviate by 0.1 or more.
hw9_4.m Predictive power (number of iterations before deviations accumulate) vs. accuracy in specifying initial conditions for the Ricker equation in chaotic regime (r=4).
5 With double precision accuracy, we cannot trust the computation yn+1=yn e4 (1 - yn) for more than 74 steps. The curve nd vs. m is approximately linear, with a slope between 4.5 and 4.9 (see plot produced by hw9_4.m). Therefore, in order to predict (calculate) the value of y500, we need between 500/5 = 100 and 500/4.7 = 106 significant digits.
Homework #10
1a Newton.m Matlab function Newton(f, df, x0, tol) calculates the root of any equation f(x)=0 with the derivative given by funtion df(x) = f '(x). To find the root of equation cos(x) = x2, using initial condition x0=1 and tolerance of 10-9, execute the command Newton( inline('x^2-cos(x)'), inline('2*x+sin(x)'), 1, 1e-9 );
NewtonSym.m Optional: same as above, but this program uses Matlab's symbolic capabilities (the Maple core) to calculate the derivative. Execute the commands x = sym('x'); Newton(x^2-cos(x), 1, 1e-9 );
1b Bisection.m Matlab function Bisection(f, a, b, tol) calculates the root of equation f(x)=0 bracketed by the interval [a,b]. To find the root of equation cos(x) = x2 between a=0 and b=π/2, with a tolerance of 10-9, execute the command Bisection( inline('x^2-cos(x)'), 0, pi/2, 1e-9 );
1c FalsePos.m Matlab function FalsePos(f, a, b, tol) calculates the root of equation f(x)=0 bracketed by the interval [a,b], using the method of false position. Execute the command FalsePos( inline('x^2-cos(x)'), 0, pi/2, 1e-9 );
1d Let's now choose the g(x) function for the iteration method: g(x) = x + C f(x) = x + C(x2 - cos x), so g'(x) = 1 + C(2 x + sin x). Stability condition is -1<g'(xroot)<1 . In fact, we know that the convergence of a difference equation is the fastest when g'(xroot) is close to zero: 1 + C(2 xroot + sin xroot) = 0, which yields C = -1 / (2 xroot + sin xroot) . Now, we know that 0<xroot<π/2 and 0<sin xroot<1. To get a rough estimate for C, let's take the values in the middle of these intervals: xroot ~ π/4, sin xroot ~ 1/2, so 2 xroot + sin xroot ~ π/2 + 1/2 ~ 2, yielding C = -1/2 (this is a crude approximation, but it's sufficient)
Iteration.m Matlab function Iteration(f, x0, tol) calculates the root of equation f(x)=0 using the method of iterations. Execute the command Iteration(inline('x-(x^2-cos(x))/2'), 1, 1e-9);
3 ConvOrder.m Matlab function ConvOrder(error) makes the plot of ln(εn+1) vs. ln(εn), using an array of error values ("error"), and calculates the order of convergence of a given method, using the first and last point on the plot. You must save this function into your Matlab directory before running the programs below
2a, 3 NewtonE.m Matlab function NewtonE(f, df, x0, tol, exact) is similar to Newton(f, df, x0, tol) (see above), but uses the exact root ("exact") to find the convergence order of the method. To find the order of convergence of the Newton method for equation x2-7=0, using initial condition x0=1 and tolerance of 10-9, execute the command NewtonE( inline('x^2-7'), inline('2*x'), 1, 1e-9, sqrt(7) );
2b, 3 BisectionE.m Function BisectionE(f, a, b, tol, exact) is similar to function Bisection(f, a, b, tol) (see above), but uses the exact root to find the convergence order of the bisection method. Execute the command BisectionE( inline('x^2-7'), 2, 3, 1e-9, sqrt(7) );
2c, 3 FalsePosE.m Same as above, but for a false position method. Execute the command FalsePosE( inline('x^2-7'), 2, 3, 1e-9, sqrt(7) );
2d, 3 Let's choose the g(x) function for the iteration method: g(x) = x + C f(x) = x + C(x2-7), so g'(x) = 1 + 2 C x. Stability condition -1 < g'(xroot) < 1 becomes -1 < 1 + 2 C xroot < 1 . Therefore, -2 < 2 C xroot < 0 , which yields -1 < C xroot < 0 . Thus, C is a negative constant satisfying C > -1/xroot = -1/sqrt(7) = -0.378. The following simple choices will work: C=-1/4, -1/5, -1/6, -1/7... The best choice is the one yielding g'(xroot)=0: C = -1 / (2 xroot) ~ -0.189 ~ -1/5
IterationE.m Function IterationE(f, x0, tol, exact) calculates the root of equation f(x)=0 and finds the order of convergence of the iteration method. Execute the command IterationE(inline('x-(x^2-7)/5'), 3, 1e-9, sqrt(7));
Homework #11
1-3 midpoint.m Function midpoint(a, b, f, N) calculates the integral of inline function f from a to b using the midpoint rule with N subintervals. Call example: midpoint(0, pi/2, inline('cos(x.^2)'), 10);
trapezoid.m Function trapezoid(a, b, f, N) calculates the integral of inline function f from a to b using the trapezoid rule with N subintervals. Call example: trapezoid(0, pi/2, inline('cos(x.^2)'), 10);
simpsons.m Function simpsons(a, b, f, N) calculates the integral of inline function f from a to b using the simpsons rule with N subintervals. Call example: simpsons(0, pi/2, inline('cos(x.^2)'), 10);
1 hw11_1.m Comparison of the three intergration methods in calculating the integral of cos(x2) from 0 to π/2.
2-3 hw11_2.m Comparison of the three intergration methods in calculating the integral of ln x from 1 to 3. Produces a plot and determines the convergence of the error.
Homework #12
1 hw12_1.m The pseudorandom congruential generator: the uniform distribution check and a check for correlations between successive numbers.
2 MC_sphere.m Function MC_sphere(N) calculates the volume of a sphere using the Monte-Carlo method with N points. Execute for i=[1000, 5000, 10000, 50000]; MC_sphere(i); end to see the results.
3 MC_ellipse.m Function MC_ellipse(N) calculates the volume of a sphere using the Monte-Carlo method with N points.
hw12_3.m Calculates the area of an ellipse given by x2/4 + y2 < 1 using the Monte-Carlo method. The exact result is A = 2π. In general, the area of an ellipse with semi-axes a and b is πab.


Lecture Notes:

Jan 19 fprintf statement: some examples
Jan 24 Making plots: an example
Feb 7 Integrating exp(-x2) using the Taylor series
Feb 9 Recursive and non-recursive exponential functions
Feb 16 Round-off errors, binary system, infinities, and NaN's
Feb 28 Solving difference equations: Fibonacci numbers
Mar 2 Review of round-off errors and loss of significance (Word doc file)
Mar 28 Bifurcation diagram for the logistic equation: yn+1 = r yn (1 - yn) (JPEG file)
Mar 30 Root-finding (Word doc file)
Apr 20 This MontyHall function performs a Monte-Carlo simulation to determine the probability of winning in the Monty Hall game, given the two different strategies. To determine the win probability for the keep strategy, execute MontyHall(1); compare this with the win probability for the switch strategy, by executing MontyHall(0)

Victor Matveev
Last modified: Feb 15, 2005
This server is running a Redhat distribution of Linux