12.010
Homework #2 Due
Tuesday, October 16, 2007
|
Questi(a) Write, compile and run a fortran program which
generates a table of Bessel functions of the first kind for integer orders
between 0 and 5, and argument between 0 and 1 in steps of 0.1. Bessel functions of the first kind
satisfy the following differential equation:
where the function y is the Bessel function of order
n. (see http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html
for information on Bessel function ). Bessel functions are rarely computed
by directly the solving the differential equation. The values in the table should be given with 5 decimal
places. The table should have
headers explaining what the columns are. Explain how you designed the program
and give an example of the output. (b) How would you change this program if 10 significant
digits were required? Fortran source code and example output should also be
supplied on (1): (25-points) |
The main decision to be made here is how
to implement the evaluation of the Bessel function. Additional considerations are the formatting of the
table but this is relatively easy in Fortran.
The http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html
gives the series expansion of the Bessel Function of
the first kind. From equation 46
we have
This formula will work for all values of the argument
x. To implement this equation we
note the sum is to infinity, which will take a long time to compute. So we need to decide on the number of
terms to include. As we look at
the expression, for |x|<1 the series is well behaved in that as l goes to
infinity, x2l+m will go to zero (i.e., the terms get smaller and
smaller). However for |x|>1, x2l+m
will go to infinity as l goes to infinity. These ever increasing values are reduced in size by the l!
and (m+l)! in the denominator. The magnitudes of l! and (m+l)! could cause
problems. (m+l)! grows rapidly and will lead to overflow of integer*4 when
there are enough terms in the series. Sample calculations of the series suggest at for
x=10.0, l+m should be 32 to obtain 10-6 accuracy in the
calculation. Factorial(32) is
2.6x1035 which 25-orders of magnitude larger than the maximum value
possible with integer*4. Therefore
although factorial in an integer, it will need to be computed as a real*8 in
this calculation.
Another aspect of this series is the alternating sign (-1)l. For arguments greater than 1, the
series will be summing and subtracting numbers that are large and this can introduce
rounding errors.
In g77 and gcc, there is an intrinsic bessel function and if
the user routine called bessel is used, a warning is printed saying that there
is a conflict between the user routine and the intrinsic routine. In g77, the intrinsic routine is used
and not the user version which is rather surprising. (Normally, user routines over rider system functions but
since this an intrinsic function conflict this rule must not apply).
Code design:
The major issue with accuracy of the series calculations and
factorial magnitude were discussed above.
The major components of the code are to
(1) Output table headers so that the user knows what is
being tabulated. In this code, the
number of orders and number of significant digits can specified and so the
formats for output are generated as strings and used to write out the results
and header.
(2) Functions to evaluate the Bessel function. In my implementation, the
besselfunction functions evaluates how many terms are needed to obtain the
accuracy set by the number of significant digits in the output (eps). The number of terms for each argument
varies because of this.
(3) Function needed to evaluate factorial. As noted above this is a real*8
function to ensure that the result does not overflow.
Solution:
The fortran code I used is linked to hw2_1_07.f.
The output of my program is given below.
G4-Computer[343]
g77 hw2_1_07.f -o hw2_1_07
G4-Computer[344]
hw2_1_07
----------------------------------
TABLE OF BESSELS ORDERS 0 TO 5
----------------------------------
--------------------------------------------------------------------------------
Arg | J00 |
J01 | J02 | J03 | J04 | J05 |
--------------------------------------------------------------------------------
0.00 | 1.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.00000 |
1.00 | 0.76520 | 0.44005 | 0.11490 | 0.01956 | 0.00248 | 0.00025 |
2.00 | 0.22389 | 0.57672 | 0.35283 | 0.12894 | 0.03400 | 0.00704 |
3.00 | -0.26005 | 0.33906 |
0.48609 | 0.30906
| 0.13203 | 0.04303 |
4.00 | -0.39715 | -0.06604 |
0.36413 | 0.43017
| 0.28113 | 0.13209 |
5.00 | -0.17760 | -0.32758 |
0.04656 | 0.36483
| 0.39123 | 0.26114 |
6.00 | 0.15065 | -0.27668 | -0.24287 |
0.11477 | 0.35764
| 0.36209 |
7.00 | 0.30008 | -0.00468 | -0.30142 |
-0.16756 | 0.15780
| 0.34790 |
8.00 | 0.17165 | 0.23464 | -0.11299 | -0.29113 |
-0.10536 | 0.18577 |
9.00 | -0.09033 | 0.24531 |
0.14485 | -0.18093 | -0.26547 | -0.05504 |
10.00 | -0.24594 | 0.04347 |
0.25463 | 0.05838
| -0.21960 | -0.23406 |
--------------------------------------------------------------------------------
How would you change this program if 10 significant
digits were required?
Several changes would be needed to output to 10-signifciant digits.
(1) The real*8 arguments would needed (already coded that way)
(2) The formats would need to change for the extra digits. Code currently does that through generating its own formats.
(3) The headings format would need to be changed to re-align the columns. Again done in code.
(4) Care is always needed that floating point constants should have d0 added to ensure that double precision values are used.
|
Question (2): (25-points).
|
|
Write a program that reads a file containing text, counts
the number of characters (letters a-z) and words in the text, and output the
text with capital letters at the beginning of each sentence. The text below is contained in the
file Q2_text.txt. we take
as a self evident foundational principle that the set of effects to be
considered as contributing to local station displacements and the conventional
models to be applied for their compensation should be guided
by rational and well considered bases, and should not be developed haphazardly
or randomly. for historical
reasons and general consistency, it might
be prudent to retain some past practices even if they are not fully
consistent with the adopted principles; but future expansions should
be determined by specified rules.
this position paper proposes such a
set of guidelines and rationales for iers Conventions updates. Hints: (1) Look at the ASCII table and check the relationship
between upper and lower case letters. Intrinsic function CHAR and ICHAR
convert between character strings and ASCII codes and visa versa. (3) Writing with the format above (instead of aMN, where MN is a
number, for example) will output only the number of characters in the string
to be output. To avoid extra
spaces, only print the number of characters needed using the 1:N feature
where N is the number of characters needed. Q2F < Q2text.txt where Q2F is the name of the program
(you can call the program any name you like). |
Again a nominally easy program but we will clearly need a few utility subroutines to perform various tasks. From the question, the main ones will be:
(1) A subroutine that will convert a string to upper case. Examination of the ASCII codes shows that the only difference between upper and lower case is that the upper case symbol for a letter is 32 less than the lower case value. (man ascii on most systems will print the ASCII table). Therefore by use of ICHAR to get the ASCII numeric value and CHAR to convert the ASCII numeric value back to a character, we can convert from lower to upper case. But we do need to be careful: Only the characters from a-z should be converted. A similar routine that converts to lower case is also useful.
(2) Since we need to do other manipulations of strings, we will probably need something that finds the length of the used portion of a string, i.e., even though a string may be declared to a certain length only some of the string is filled with non-blank characters (in class we discussed that fortran pads strings with blanks if a shorter length string is assigned to a longer string. F90 has an intrinsic that does this all lentrim but this routine
(3) Routines are needed that counts characters (a-z) (easily done) and counts works (a little trickier)
(4) Routine to capitalize starts of sentences. Checks for punctuation and then waits for next (a-z) character which might be on the next line from the file.
For items 2 and 3 above, the routines will work by finding the positions in strings of blank and non-blank characters and using the fortran string(n:m) feature to extract those parts of the string.
Code design
á Ask the user for the string that contains the file name to be read. By pre-loading the string with a default file name and then read a free format read (*) so that /return will keep the default value.
á Check for IOSTAT errors on opening the file status 'old'
á Loop over reading lines from file until an error reading file (normally error will be -1 meaning end-of-file reached).
á Count characters and words
á Capitalize as needed in output updated line
á At end of program output final counts.
Routine tocaps( in_string, out_string)
Get the minimum length of string passed (so that we donÕt exceed bounds)
Loop over each character in string
If the input character is between ÔaÕ and ÔzÕ
Subtract 32 from ASCII code and convert back to character
Endif
End loop
Function lenline(inline)
Returns the length of the used portion (non-blank character) of a string
Get the full length of the string
Start at the end of the string working backwards to find a non-blank character. Make sure we do not go past the start of the string
Return the number of characters to last non-blank part of the string
The Fortran code I used is linked to hw2_2_07.f
G4-Computer[345] hw2_2_07
12.010 HOMEWORK 2 QUES 2:
Read and count file
Enter file name (</cr>
for default Q2_text.txt) /
Formatted version of
Q2_text.txt
We take as a self evident
foundational principle that the set of effects
to be considered as
contributing to local station displacements and the
conventional models to be
applied for their compensation should be
guided by rational and well
considered bases, and should not be developed
haphazardly or
randomly. For historical reasons
and general consistency,
it might be prudent to
retain some past practices even if they are not
fully consistent with the
adopted principles; but future expansions
should be determined by
specified rules. This position
paper proposes
such a set of guidelines and
rationales for iers conventions updates.
File Q2_text.txt has 533 characters and 99 words
|
Question (3): (50-points) Write a Fortran program that will compute the motion of a baseball in 2-dimensions. The program should generate the trajectory of the ball and compute the change in height of the ball from where it is initially released and home plate 18.44 m away and the velocity at home plate. The dynamical forces acting on the baseball are:
where Ar
is the cross-sectional area of the ball, Cd
is the drag coefficient which depends on the velocity, v, of the ball, r is the density of air.
The drag coefficient depends on velocity with vd =35 m/s and D=5m/s.
The Magnus force, Fm
is the lift force due to the rotation of the ball w. S is a
scaling coefficient for the Magnus force with m being the mass of the baseball. Although S probably depends on velocity, in this
problem we will assume that it is constant. (see discussion in: http://farside.ph.utexas.edu/teaching/329/lectures/node43.html
) In your problem, the rotation vector will be perpendicular to the velocity vector (i.e., horizontal) and hence the Magnus force will only lift or depress the ball and not move it out of its plane of motion. As a test of your program use the following constants to compute the trajectory of the ball: Assume the following values Baseball circumference 23 cm Baseball mass (m) 150g r = 1.226 kg/m3 g = 9.8 m/s2 Initial release angle +1 deg Initial velocity 85 mph Spin rate 1800 rpm (counter clockwise so ball will lift). The height at home plate should be calculated to an accuracy of 1mm. Hint: Be very careful with units (there are all mixed above). Your answer to this question should include: (a) The algorithms used and the design of your program (b) The Fortran program source code (I will compile and run your programs). (c)
The results from the test case above that will include the
change in height from release to home plate (18.44 m) and the velocity at
home plate. |
With these equations, we can compute the forces acting on the plane (gravity, Magnus force and drag) and using the mass of the ball, compute accelerations. By integrating the accelerations we generate velocity changes and by integrating the velocity changes we get position changes. Since the forces depend on velocity, we need to integrate in small steps. Precisely how small we discuss below. This problem is basically an integration problem. It can also be posed as the solution to a second-order differential equation and in later homework (and languages) we will how to solve it this way.
The first thing to do in a problem of this type is go to the library and read about numerical integration. Notice that there are two parts to this problem: (1) you need to do the integration and (2) the program has to determine the accuracy of the integration. In particular, if the user requests low accuracy integration, then the program should run quicker. In reading about numerical integration both methods for integration and the way of judging accuracy must be determined. You need to know the techniques well enough to implement in a program.
As discussed in the solution to Homework 1, the basic equation being integrated here the acceleration of the body. The tricky part of the problem is that the acceleration depends on the velocity of the body. We are solving a pair of non-linear ordinary differential equations in this problem. We have:
where x is a vector of the x and z coordinates of the object, and the derivatives are with respect to time. The k term here is lift and drag. If k is zero, then the equation reduces to the simple one that can be analytically solved.
I consulted Numerical Recipes (for the discussion of numerical integration, not the programs) and Abramowitch and Stegum. The latter book gives the Runge-Kutta formula for numerically integrating an equation of the form yÕÕ(x) = f(x,y,yÕ) where in our case x is time, y is position, yÕ is velocity and yÕÕ is acceleration. (You can also look at Lecture 24 of the OCW version of 12.010).
The equations above give the position (y) and velocity (yÕ) at time step n+1 in the numerical integration based on the values at time step n. The length of the time step is h. The O(h^5) statement means that the error in the integration decreases as step size to the 5th power. The simple approach to the numerical integration in EulerÕs method in which the acceleration (function f) is computed for the position and velocity at the start of the step and is assumed constant for the distance traveled and velocity change during the step. The above equation is an approach (one of many different ones available) in which the constant acceleration assumption is replaced with accelerations computed at different points in the time step. A very simple application of this approach is consider where gravity will be computed. Since the velocity of the body is known at the start of the step, an approach height at the end of the step can be computed (assuming constant velocity) and constant value of gravity used could be computed at the mid-point of the height at the beginning of the step and the approximate position at the end of the step. Similarly, rather than computing drag using the initial velocity only, the change in velocity can be estimated and the drag computed using the average velocity of the drag during the step. The Runge-Kutta equations are derived using these ideas.
The Runge-Kutta integration method allows the trajectory to be computed based in an initial velocity. Our problem is two-dimensional and in f90 we could implement the equations above uses two element arrays for y and yÕ. In f77, we could compute each component separately, or we could use a technique commonly used in 2D problems. We can make y and yÕ complex data types where the real component is the x-coordinate and the imaginary component is the z coordinate.
One complication in the integration here is that we are not integrating for a finite length of time but rather until the z-coordinate reaches zero. This is a bit tricky for the last integration step where the time interval until the objects is not likely to be the time step being used for the integration. The implementation in the solution is to integrate until the object goes below the surface. Based on how far back in time in that last step, the last step of the integration is done with a smaller step size. The last step is repeated with smaller step sizes until the object is within the tolerance of the surface for the accuracy of the calculation.
The basic flow of the program is to:
-Read the inputs fro the user. These reads are ordered such that those quantities most likely to be changed by the user. These values have defaults and so the user only needs to use /<cr> to adopt the defaults. - Compute the initial velocity and derivative based on the analytic solution
- Trial integrations with decreasing step sizes to determine the step size needed to hit the ground within a specific accuracy.
- Note the code that tests for the ball arriving at home plate. The integration step size is reduced so that we converge onto home plate
- Once the step size is set, run with the final step size and output results. (This later step is a little inefficient since it repeats an earlier calculation but with output generated this time.
Out put the final height change to home plate and the velocity at home plate.
The program in implemented is given in hw2_3_07.f. It also uses an include file BBall.h that contains parameters for gravity values, air density etc and a common block that contains the Òfixed:Ó variables in the program (object mass, etc.)
Example output for test case:
G4-Computer[347] hw2_3_07
12.010 Program to solve
Baseball Trajectory problem
Given initial velocity and
spin, program computes the trajectory of the ball and the height change
to reach homeplate
Program Parameter Input.
[Defaults are printed use /<cr> to accept default
Pitch Speed [85.00 mph],
Pitch Angle [1.0deg], Spin Rate [1800.0 rpm]
Pitch velocity: 38.00 m/s, Angle 1.0 deg, Spin Rate 188.50
rads/s
MASS 0.1500 kg, Area 0.004 m^2
GRAVITY 9.80000 m/s**2 RHO AIR 1.22500 kg/m**3
Initial Velocity 37.998
m/s;
85.000 mph
Step size 0.100
sec
Step 0.05
-0.566211262 -0.556799502
9.4117603
Step 0.025
-0.538202987 -0.538709154
0.506167751
Baseball Trajectory
T Time (s) X_pos (m) Z_pos (m) X_vel (m/s) Z_vel (m/s) Step (s)
S 0.00000
0.00000 0.00000
37.9926 0.6632 0.0250000
F 0.02500
0.94715 0.01439
37.7794 0.4881 0.0250000
F 0.05000
1.88897 0.02441
37.5666 0.3137 0.0250000
F 0.07500
2.82548 0.03008
37.3540 0.1398 0.0250000
F 0.10000
3.75668 0.03140
37.1418 -0.0336 0.0250000
F 0.12500
4.68257 0.02840
36.9299 -0.2063 0.0250000
F 0.15000
5.60318 0.02109
36.7184 -0.3785 0.0250000
F 0.17500
6.51850 0.00948 36.5073 -0.5501 0.0250000
F 0.20000
7.42855 -0.00641
36.2967 -0.7210 0.0250000
F 0.22500
8.33334 -0.02656
36.0865 -0.8914 0.0250000
F 0.25000
9.23288 -0.05097
35.8768 -1.0611 0.0250000
F 0.27500
10.12718 -0.07961
35.6676 -1.2302 0.0250000
F 0.30000
11.01626 -0.11248
35.4590 -1.3987 0.0250000
F 0.32500
11.90013 -0.14954
35.2509
-1.5666 0.0250000
F 0.35000
12.77881 -0.19080
35.0435 -1.7338 0.0250000
F 0.37500
13.65232 -0.23623
34.8368 -1.9003 0.0250000
F 0.40000
14.52066 -0.28581
34.6307 -2.0662 0.0250000
F 0.42500
15.38386 -0.33953
34.4254 -2.2315 0.0250000
F 0.45000
16.24193 -0.39738
34.2208 -2.3961 0.0250000
F 0.47500
17.09491 -0.45933
34.0171 -2.5600 0.0250000
F 0.50000
17.94279 -0.52538
33.8141 -2.7233 0.0250000
F 0.50981
18.27417 -0.55241
33.7347 -2.7872 0.0098116
F 0.51309
18.38468 -0.56158
33.7082 -2.8085 0.0032771
F 0.51418
18.42156 -0.56465
33.6994 -2.8156 0.0010941
F 0.51455
18.43385 -0.56568
33.6965 -2.8180 0.0003649
F 0.51467
18.43795 -0.56603
33.6955 -2.8188 0.0001216
F 0.51471
18.43932 -0.56614 33.6952
-2.8190 0.0000406
F 0.51472
18.43977 -0.56618
33.6950 -2.8191 0.0000135
F 0.51473
18.43992 -0.56619
33.6950 -2.8192 0.0000045
A 0.51473 18.44008 -0.56620
33.6950 -2.8192 0.0000015
Final Height change -.5662
m; Velocity 75.6 mph Angle -4.78 deg