Final Exam  Matlab Part  
1  stephensen.m 
Matlab function stephensen(f, x0, tol) calculates the root of equation
f(x)=0 , using initial condition x_{0} .
execute stephensen( inline('exp(x)x'), 1, 1e8 );

2  dice.m 
Matlab function dice(diceSum, nTrials) calculates the average number of twodice rolls
required to get a sum of diceSum , using a MonteCarlo 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)

Homework #3  
1  hw3_1.m  Limit of (e^{x}1) / x for x > 0 
2  hw3_2.m  Limit of (root(x^{2}+1)  1) / x^{2} for x > 0 
3  IntCosX2.m  Matlab function integrating cos(x^{2}) 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 x^{2} ((x1)^{1}(x+1)^{1}) and g(x)=x^{2} (x^{2}  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 roundoff 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=1x 
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 10^{15} ( a_{7} = x^{8} y^{7} 15! / (7! 8!) = 1.04 * 10^{16} ), 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 10^{15}. When n=14, the k=7 term is a_{7} = 10^{7} 9^{7} 14! / (7! 7!) = 1.64 * 10^{17}  
3  mySin.m  Matlab function calculating sin(x) by summing the Taylor series 
hw4_3.m  Table comparing mySin(x) to builtin 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 3540 radian. This is because the terms in the Taylor series are becoming very large as x increases (larger than 10^{15}), leading to roundoff 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 roundoff error in representing πn becomes comparable to
x . The roundoff 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=10^{m}+1 , and m=1:17 , for the builtin 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 1516.


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 builtin 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 1e15 to 1e16 , since we know
that such a small number will get completely rounded off when added to 1 . Now, for x 
Homework #6  
1  hw6_1.m  Solution is y_{n}=(2/3)^{2n}

2  hw6_2.m  Solution yields r_{1}=1/2, r_{2}=3/2, c_{1}=2, c_{2}=0 ,
so y_{n}=2/2^{n}=1/2^{n1}

3  hw6_3.m  Constant solution of nonhomogeneous equation is C = 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 roundoff errors (at least untill y_{n} become larger than 2^{53}) 
Homework #7 
Homework #8  
1a  hw8_1a0.m 
y_{n+1} = 64 y_{n} / (1 + y_{n})^{3} .
Equilibria: y_{eq}=0 and y_{eq}=3 .
Stability: f '(y) = 64 (1  2 y) / (1 + y)^{4}
(monotonically unstable),
y_{eq}=3: f '(3) = 5/4 (oscillatorily unstable);

hw8_1a3.m  
1b  hw8_1b.m 
y_{n+1} = y_{n} ln(y_{n}) .
Equilibria: y_{eq}=0 and y_{eq} = e = 2.718... .
Stability of nonzero equilibrium: f '(y) = ln(y) + 1 > f '(e) = 1 + 1 = 2
(monotonically unstable)

2 
y_{n+1} = y_{n} e^{r (1  yn)} .
Equilibria: y_{eq}=0 and y_{eq} = 1 .
Stability: f '(y) = (1  r y) e^{r (1  yn)} Zero equilibrium: f '(0) = e^{r} : stable if r < 0 ;
Nonzero 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 y_{0}= 
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 y_{n+1}=y_{n} e^{r (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 y_{n+1}=y_{n} e^{4 (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 y_{500} , 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) = x^{2} , using initial condition
x_{0}=1 and tolerance of 10^{9} , execute the command
Newton( inline('x^2cos(x)'), inline('2*x+sin(x)'), 1, 1e9 );

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^2cos(x), 1, 1e9 );


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) = x^{2} between a=0 and b=π/2 
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^2cos(x)'), 0, pi/2, 1e9 );

1d 
Let's now choose the g(x) function for the iteration method: g(x) = x + C f(x) = x + C(x^{2}  cos x) ,
so g'(x) = 1 + C(2 x + sin x) . Stability condition is 1<g'(x_{root})<1 .
In fact, we know that the convergence of a difference
equation is the fastest when g'(x_{root}) is close to zero: 1 + C(2 x_{root} + sin x_{root}) = 0 ,
which yields C = 1 / (2 x_{root} + sin x_{root}) . Now, we know that 0<x_{root}<π/2
and 0<sin x_{root}<1 . To get a rough estimate for C, let's take the values in the
middle of these intervals: x_{root} ~ π/4, sin x_{root} ~ 1/2, so 2 x_{root} + sin x_{root} ~ π/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^2cos(x))/2'), 1, 1e9);


3  ConvOrder.m 
Matlab function ConvOrder(error) 
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 x^{2}7=0 , using initial condition
x_{0}=1 and tolerance of 10^{9} , execute the command
NewtonE( inline('x^27'), inline('2*x'), 1, 1e9, 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^27'), 2, 3, 1e9, sqrt(7) );

2c, 3  FalsePosE.m 
Same as above, but for a false position method. Execute the command
FalsePosE( inline('x^27'), 2, 3, 1e9, sqrt(7) );

2d, 3 
Let's choose the g(x) function for the iteration method: g(x) = x + C f(x) = x + C(x^{2}7) ,
so g'(x) = 1 + 2 C x . Stability condition 1 < g'(x_{root}) < 1 becomes
1 < 1 + 2 C x_{root} < 1 . Therefore, 2 < 2 C x_{root} < 0 , which yields
1 < C x_{root} < 0 . Thus, C is a negative constant satisfying C > 1/x_{root} = 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'(x_{root}) 

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^27)/5'), 3, 1e9, sqrt(7));

Homework #11  
13  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(x^{2})
from 0 to π/2.

23  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 MonteCarlo 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 MonteCarlo method
with N points.

hw12_3.m 
Calculates the area of an ellipse given by x^{2}/4 + y^{2} < 1
using the MonteCarlo method. The exact result is A = 2π . In general,
the area of an ellipse with semiaxes a and b is πab .

Jan 19  fprintf statement: some examples 
Jan 24  Making plots: an example 
Feb 7  Integrating exp(x^{2}) using the Taylor series 
Feb 9  Recursive and nonrecursive exponential functions 
Feb 16  Roundoff errors, binary system, infinities, and NaN's 
Feb 28  Solving difference equations: Fibonacci numbers 
Mar 2  Review of roundoff errors and loss of significance (Word doc file) 
Mar 28  Bifurcation diagram for the logistic equation: y_{n+1} = r y_{n} (1  y_{n}) (JPEG file) 
Mar 30  Rootfinding (Word doc file) 
Apr 20  This MontyHall function performs a MonteCarlo 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) 