## Programming Exercises

The following exercises are meant to be answered by either a script or a function Mfile (you decide). The descriptions are complete but the nature of the input/output and display options is left to you. Also, some problems are much simpler than others. Suggestion: start simple and add complexity.

[Back to main page] [Back to exercises page]
```1. Compute the value of pi using the following series

How many terms are needed to obtain an accuracy of 1e-12?  How accurate is
the sum of 100 terms of this series?

2. The Fibonacci numbers are comuted according to the following relation:

Fn = Fn-1 + Fn-2

with F0 = F1 = 1.

a. Compute the first 10 Fibonacci numbers.
b. For the first 50 Fibonacci numbers, compute the ratio

Fn / Fn-1

It is claimed that this ratio approaches the value of the golden mean
( (1 + sqrt(5))/2 ).  What do your results show?

3. The Legendre polynomials (Pn(x)) are defined by the following recurrance relation

(n+1) Pn+1(x) - (2n+1)x Pn(x) + n Pn-1(x) = 0

with P0(x) = 1, P1(x) = x and P2(x) = (3x2 - 1)/2. Compute the next three
Legendre polynomials and plot all 6 over the interval [-1,1].

4. The present value of an annuity (a yearly sum of money) may be computed from the
formula

P = (A/i) [(1 + i)n - 1]/(1 + i)n

where A is the annuity (in \$/year), i is the nominal yearly interest rate (in
decimal form), n is the number of years over which the annuity is paid and P is
the present value (\$).

Example computation: If i = 0.15 (15%), A = \$100/year and n = 10 years then
P = \$501.88.

If you won the \$1,000,000 State Lottery and the Lottery offered you the choice of
\$500,000 today or \$50,000/year for 20 years, which would you take?  You can assume
an inflation (interest) rate of 5%.

5. I found the following algorithm on a web page* for computing pi:

1. Set a = 1, b = 1/sqrt(2), t = 1/4 and x = 1

2. Repeat the following commands until the difference between a and b
is within some desired accuracy:

y = a
a = (a + b)/2
b = sqrt(b*y)
t = t - x*(y - a)^2
x = 2*x

3. From the resulting values of a, b and t, an estimate of pi is

Pi_est = ((a + b)^2)/(4*t)

How many repeats are needed to estimate pi to an accuracy of 1e-8? 1e-12?
Compare the performance of this algorithm to the result from Exercise 1, above.

*/Gauss_L.html

6. Write a script that asks for an integer (n) and then computes the following based
on the value of the integer:

While the value of n is greater than 1, replace the integer with
half of its value (n/2) if the integer is even.  Otherwise, replace
the integer with three times its value, plus 1 (3*n + 1).

Make provision to count the number of values in (or the length of) the sequence
that results.

Example calculation: If n = 10, the sequence of integers is 5, 16, 8, 4, 2, 1 and
so the length is 6.

Make a plot of the length of the sequence that occurs as a function of the
integers from 2 to 30.  For example, when n = 10, the length is 6 while for
n = 15, the length is 17.  Is there any pattern?  Try larger numbers to see
if any pattern occurs.  Is there any integer for which the sequence does
not terminate?

7. Write a script/function that converts a Roman numeral to its decimal equivalent.
There are two distinct situations that you might design your program for:

a. The "old" style where the order of the symbols does not matter.  In this
case, IX and XI both mean 10 + 1 or 11.  You should be able to handle the
following conversion table:

Roman          Decimal
I                1
V                5
X               10
L               50
C              100
D              500
M             1000

b. The "new" style where the order of the symbols does matter.  For example,
IX is 9 (10 - 1), XC is 90 (100 - 10).  The conversion table given above
still holds and you may assume for this case that the only instances of
"order" you will encounter are

IV (4), IX (9), XL (40), XC (90), CD (400) and CM (900)

The function input will be useful here.  The format

>> str = input('Roman numeral: ','s')

will provide a way to get the Roman number into your program as a string.  It
would be a good idea to try case a. first.

8. Write a function that will do the inverse of the previous problem - convert a
decimal number into a Roman number.

9. Compute and plot the path(s) of a set of random walkers which are confined by
a pair of barriers at +B units and -B units from the origin (where the walkers
all start from).

A random walk is computed by repeatedly performing the calculation

xj+1 = xj + s

where s is a number drawn from the standard normal distribution (randn in
MATLAB).  For example, a walk of N steps would be handled by the code fragment

x(1) = 0;
for j = 1:N
x(j+1) = x(j) + randn(1,1);
end

There are three possible ways that the walls can "act":

a. Reflecting - In this case, when the new position is outside the walls, the
walker is "bounced" back by the amount that it exceeded the barrier.
That is,

when xj+1 > B,
xj+1 = B - |B - xj+1|

when xj+1 < (-B),
xj+1 = (-B) + |(-B) - xj+1|

If you plot the paths, you should not see any positions that are beyond |B|
units from the origin.

b. Absorbing - In this case, if a walker hits or exceeds the wall positions, it
"dies" or is absorbed and the walk ends.  For this case, it is of interest
to determine the mean lifetime of a walker (i.e., the mean and distribution
of the number of steps the "average" walker will take before being absorbed).

c. Partially absorbing - This case is a combination of the previous two cases.
When a walker encounters a wall, "a coin is flipped" to see if the walker
relfects or is absorbed.  Assuming a probability p (0 < p < 1) for
reflection, the pseudo-code fragment that follows uses the MATLAB uniform
random-number generator to make the reflect/absorb decision:

if rand < p
reflect
else
absorb
end

What do you do with all the walks that you generate?  Compute statistics, of course.

What is the average position of the walkers as a function of time?
What is the standard deviation of the position of the walkers as a function
of time?
Does the absorbing or reflecting character influence these summaries?
For the absorbing/partial-reflection case, a plot of the number of surviving
walkers as a function of step numbers is a very interesting thing.

is useful, informative and interesting, particularly if graphically displayed.

10. The properties of saturated steam are tabulated in many books but that is not
a useful form for computer use.  The following expressions were provided by
Williamson (Chemical Engineering, May 15, 1972, p. 128):

Saturation temperature, degrees F

Tsat = 8576.65/(15.47538 - ln(Psat)) - 459.216 - 0.023719Psat +
(0.84219e-4)Psat2 - (0.70854e-7)Psat3

Liquid specific volume, ft3/#

Vliq = 0.01655 + (0.150326e-4)Psat - (0.40488e-7)Psat2 +
(0.665584e-10)Psat3 - (0.4053e-13)Psat4

Vapor specific volume, ft3/#

Vvap = 430.8419/(Psat + 1.66) + 0.2031 - (0.000258)Psat

Liquid specific enthalpy, BTU/#

Hliq = 6473.878/(14.01875 - ln(Psat)) - 391.6036 + (0.022915)Psat

Vapor specific enthalpy, BTU/#

Hvap = 1142.342 + (0.76833)Psat - (0.004194)Psat2 + (0.11642e-4)Psat3 -
(0.157e-7)Psat4 + (0.8086e-11)Psat5

In all cases, the pressure is in psia and the equations are valid for pressures
between 20 and 600 psia.  Maximum error in any computation is below 1%.

a. Create a function that provides these steam properties for a given vector of
pressures.  To save typing, you should be able to copy directly from the
browser window and paste into the MATLAB editor.  Some further editting
will be required, though.

b. Write a script that allows a user to get different steam properties as
desired.  The script should run until the user decides to quit.

c. Extend your script to allow the user to specify the units set (e.g., SI,
English, CGS) before requesting numerical values for the properties.

d. Extend the function from a. to provide the heat of vaporization of water
as a property.

e. What if the user wants to specify the temperature rather than the pressure?

11. Write a function which computes the cumulative product of the elements in a
vector.  The cumulative product of the jth element of the vector x, xj, is defined by

pj = (x1)(x2) ... (xj)

for j = 1:length of the vector x.  Create 2 different versions of this function:

a. One that uses two for-loops to explicitly carry out the calculations,
element by element.  An "inner" loop should accumulate the product and an
"outer" loop should more through the elements of the vector p.

b. One that uses the built-in function prod to replace the inner loop.

In each case, you can check your results with the built-in function, cumprod.

12. Follow the directions for Exercise 11 but create a function that computes the
cumulative sum of the elements of a vector.  The elements of the cumulative sum
vector are defined by

sj = x1 + x2 + ... + xj

for j = 1: length of the vector x.

The built-in functions sum and cumsum should be used instead of prod and cumprod,
respectively.

13. Create a function that generates random arrays of integers beween a and b,
inclusive.  Use the definition line

function A = randint(a,b,M,N)

where a and b define the (inclusive) range and M and N define the size
of the output array (rows and columns, respectively).

a. Test your function with the following piece of code:

x = randint(10,17,100000,1);
hist(x,10:17)

The resulting histogram should be nearly flat from 10 to 17.  Pay particular
attention to the endpoints of the distribution.

b. Test your function with the following pieces of code:

x = randint(-30,5,100000,1);
hist(x,-30:5)

x = randint(-45,-35,100000,1);
hist(x,-45:-35)

x = randint(7,-2,100000,1);
hist(x,-2:7)

Consider that each of these uses is valid.

c. Modify your code to allow for missing/default inputs.  Use the behavior of
rand to help design your code.  You will need to use nargin in your
routine.  For example, the function should be able to return

- a 5-by-5 array of integers between 1 and 20:  e.g., A = randint(1,20,5)
- a single random integer between 10 and 50:    e.g., A = randint(10,50)
- an empty array if not enough inputs are provided.

```

[Back to main page] [Back to exercises page]
Comments? Contact Jim Maneval at maneval@bucknell.edu