GEOMASS BRIEF ABSTRACT Determine the gravitational force at a point from a nearby mass and determine the underlying mass distribution given force information. GENERAL INFORMATION FileName: GEOMASS Full title: Determining the mass distribution of a bar from observed forces using physical principles. Last Revision Date: 29 May 1996. Developer: Brian J. Winkel, Department of Mathematical Sciences, United States Military Academy, West Point NY 10996 USA. Phone: 914-938-3200. Email: ab3646@usma2.usma.edu. FAX: 914-938-2409. Contact: Brian J. Winkel, Department of Mathematical Sciences, United States Military Academy, West Point NY 10996 USA. Phone: 914-938-3200. Email: ab3646@usma2.usma.edu. FAX: 914-938-2409. Aaron D. Klebanoff, Department of Mathematics, Rose-Hulman Institute of Technology, Terre Haute IN 47803 USA. Phone: 812-877-8151. Email: Aaron.Klebanoff@Rose-Hulman.Edu. FAX: 812-877-3198. Support: The production of this material is supported by the National Science Foundation under Division of Undergraduate Education grant DUE-9352849: Development Site for Complex, Technology-Based Problems in Calculus with Applications in Science and Engineering and the Arvin Foundation of Columbus IN. STATEMENT OF PROBLEM Introduction "The general problem of geological prospecting is to determine the location, shape, and constitution of subterranean bodies from measurements at the earth's surface. We will treat a much simplified one-dimensional model. Namely, we suppose that on a horizontal line (of length 1 m) measurements are made of the vertical component of gravitational force due to a variable distribution of mass along a parallel line one unit (m) below the surface." [1, p. 9] From Newton's law of gravitation it is known that there is a force of attraction between two masses, say m1 and m2 (in kg), at a distance of r m apart equal to F = G (m1 m2)/r^2 and the direction of the force is along the line connecting the masses. G is the gravitational constant in N m^2/kg^2 or m^3/(kg s^2). Here G = 6.67 x 10^(-11) m^3/(kg s^2). Forward Problem - Problem 1 We have the situation (see sketch below) in which there is a bar (1 m long) of mass density x(t) = 100 + 90 t^3 kg/m running along the x-axis, beginning at the origin (0, 0) and going to (1, 0), i.e., t goes from 0 to 1. We consider another unit point mass (1 kg) at a point along the line y = 1, say (s, 1). Distances in meters. Determine the vertical component of the gravitational force between the bar and the unit point mass (1 kg) particle at each point (s, 1) (0 <= s <= 1). Inverse Problem Approach - Problems 2 and 3 Problem 2 The task here is to determine the underlying mass density, x[t] kg/m^3, of the unit long mass at distance t meters from left hand edge - see diagram above. Observations about the vertical component of the gravitational force on a unit point mass of 1 kg which is resting 1 m above the bar where the particle is distance s m from left hand edge are known. Using an understanding of gravitational force, components of the vector of force, and addition of elements of force with the concept of integration determine the parameters of the model x(t) = a + b t which fit the data best. This then would be a reasonable conjecture as to what the underlying mass density really is from the observed surface data. Suppose we have the following experimental data for the magnitude of the vertical component of the force of attraction between the bar and the unit point mass at (s, 1). The first coordinate is s and the second coordinate is the magnitude of the vertical component of the force. data = {{0, 5.8*10^-11}, {0.2, 6.8*10^-11}, {0.4, 7.4*10^-11}, {0.6, 7.6*10^-11}, {0.8, 7.1*10^-11}}; Problem 3 Suppose we have the following experimental data for the magnitude of the vertical component of the force of attraction between the bar and the unit point mass at (s, 1). The first coordinate is s and the second coordinate is the magnitude of the vertical component of the force. data = {{0, 5.1*10^-11}, {0.2, 5.9*10^-11}, {0.4, 6.4*10^-11}, {0.6, 6.3*10^-11}, {0.8, 5.8*10^-11}}; Using an understanding of gravitational force, components of the vector of force, and addition of elements of force with the concept of integration determine the parameters of the model x(t) = a + b t + c t^2 which fit the data best. This then would be a reasonable conjecture as to what the underlying mass density really is from the observed surface data. KEYWORDS Newton's law of gravitational force of attraction between masses, vectors, vector projection, integration, function of integral, least squares, parameter estimation, inverse problem. TEACHER NOTES ISSUES RELATED TO THE PROBLEM Prerequisites Concurrent or previous exposure to Newton's law of gravity, notion of force, force vector, projection of force vector, integration as a process of accumulation, least squares optimization for parameter estimation. Time allotment - time management With some guidance this problem can be given as a homework problem with about two days to give students chance to clarify their approach. Introducing it in class should take 5 minutes the first day with another 10 minutes the second day for clarifications. Expectations We expect students to see the difference between the forward and inverse problems. Expect students to have difficulty with finding and adding up the vertical components of the forces. Future payoffs Students will have seen an inverse problem and just how difficult it must be to solve the general inverse problem by doing some suggested modiffied inverse problems. Extensions In the forward problem: One could offer several experimental values for vy[s], say at s = 0, .2, .4, .6, .8, and 1.0. Then after offering up several models of what the underlying mass densities x[t] could be one could ask which model for x[t] is most plausible given the data. One could again pose several experimental values for vy[s], say at s = 0, .2, .4, .6, .8, and 1.0 and offer up one model with an underlying set of parameters, say a quadratic x[t] = a t^2 + b t + c, and ask what are the best fitting parameters? In the inverse problem: Certainly the students could by now understand the very difficult nature of the inverse problem, namely, given the observations, what form and exactness could we offer for the underlying mass density distribution? Moreover, this approach does give some insight into the nature of the inverse problem in this case, while actually honing skills from the forward problem. References and Sources 1. Groetsch, Charles W. 1993. Inverse Problems in the Mathematical Sciences. Weisbaden GERMANY: Vieweg Publishers. pp. 9-10. 2. Halliday, David, Robert Resnick, and Jearl Walker. 1993. Fundamentals of Physics, Fourth Edition. New York: John Wiley & Sons. Chapter 15: Gravitation. POSSIBLE SOLUTION(S) Problem 1 G is the gravitational constant in N m^2/kg^2 or m^3/(kg s^2). G = 6.67 10^(-11); Let x[t] be the mass density kg/m^3 of the one meter long mass at distance t meters from left hand edge. x[t_] := 100 + 90t^3 Then vy[s] is the vertical component of the gravitational force on a particle of constant mass density, say 1, which is resting 1 m above the bar where s is distance in m from left hand edge. We obtain an integral by first determining the little element of force from (s, 1) to (t, 1) as DF = G 1 x(t) dt/(Sqrt[(s-t)^2 + 1])^2 . We then take the vertical component of this force, obtained by multiplying this magnitude DF by the cosine of the angle between the vertical and the line of direction for this force from (s, 1) to (t, 1), and integrate this over the entire line from (0,0) to (1, 0). vy[s_] = G Integrate[((s-t)^2 + 1)^(-3/2) x[t],{t,0,1}] 2 4 -11 -10 (18 - 10 s + 9 s - 9 s ) 6.67 10 (----------------------------- + 2 Sqrt[1 + s ] 2 3 4 10 (37 - 55 s + 9 s + 9 s - 9 s ) ----------------------------------- + 2 Sqrt[2 - 2 s + s ] 270 s ArcSinh[1 - s] + 270 s ArcSinh[s]) We plot a sketch of vy(s) noting its necessary asymmetry due to the non-uniform distribution of mass in our bar. Plot[vy[s],{s,0,1}, AxesLabel->{"Distance-s","MagVertForce"}] -Graphics- Problem 2 We have the following experimental data for the magnitude of the vertical component of the force of attraction between the bar and the unit point mass at (s, 1). The first coordinate is s and the second coordinate is the magnitude of the vertical component of the force. data = {{0, 5.8*10^-11}, {0.2, 6.8*10^-11}, {0.4, 7.4*10^-11}, {0.6, 7.6*10^-11}, {0.8, 7.1*10^-11}}; We plot the data: ysp = ListPlot[data,PlotStyle->PointSize[.03]] -Graphics- We now attempt a model function x(t) to see if we can "determine" the parameters which fit the z(s) function best. x[t_] = a + b t; We compute the vertical component of the force for this function x(t). yy[s_] = G Integrate[((s-t)^2 + 1)^(-3/2) x[t],{t,0,1}] 2 2 -11 a - b - a s + b s - b s b + a s + b s 6.67 10 (------------------------ + --------------) 2 2 Sqrt[2 - 2 s + s ] Sqrt[1 + s ] We now determine the least square function of the parameters a and b in our function x(t) = a + b t. We seek to minimize ss(a, b) . ss[a_,b_] = Sum[(yy[data[[i,1]]] - data[[i,2]])^2, {i,1,Length[data]}]//N; We use Mathematica's FindMinimum command to determine the values of a and b which fit the z(s) data best. sol = FindMinimum[ss[a,b],{a,1},{b,1}] -25 {3.66547 10 , {a -> 1.01831, b -> 0.503131}} We create our best fitting model function. xs[t_] = x[t]/.sol[[2]] 1.01831 + 0.503131 t So how good is our model for x(t), i.e., xs(t)? We put it into our z(s) function and plot the theoretical z's against the observed ones. z[s_] = G Integrate[((s-t)^2 + 1)^(-3/2) xs[t],{t,0,1}] 2 -11 -0.503131 - 1.01831 s - 0.503131 s 6.67 10 (-(-----------------------------------) + 2 Sqrt[1 + s ] 2 0.515182 - 0.515182 s - 0.503131 s -----------------------------------) 2 Sqrt[2 - 2 s + s ] zplot = Plot[z[s],{s,0,1}] -Graphics- Not all that bad!!!!! Show[ysp,zplot] -Graphics- Problem 3 Suppose we have the following experimental data for the magnitude of the vertical component of the force of attraction between the bar and the unit point mass at (s, 1). The first coordinate is s and the second cooredinate is the magnitude of the vertical component of the force. data = {{0, 5.1*10^-11}, {0.2, 5.9*10^-11}, {0.4, 6.4*10^-11}, {0.6, 6.3*10^-11}, {0.8, 5.8*10^-11}}; Then ys[s] is the vertical component of the gravitational force on a particle resting 1 m above the bar where s is distance in m from left hand edge. ys[s_] := G Integrate[((s-t)^2 + 1)^(-3/2) x[t],{t,0,1}] We plot the data: ysp = ListPlot[data,PlotStyle->PointSize[.03]] -Graphics- We now attempt a model function x(t) to see if we can "determine" the parameters which fit the ys(s) function best. x[t_] = a + b t + c t^2 2 a + b t + c t We compute the vertical component for this function x(t). zz[s_] = G Integrate[((s-t)^2 + 1)^(-3/2) x[t],{t,0,1}] -11 6.67 10 ( 2 2 3 a - b - c - a s + b s - c s - b s + c s - c s ------------------------------------------------\ 2 Sqrt[2 - 2 s + s ] 2 3 b + a s + c s + b s + c s + --------------------------- + 2 Sqrt[1 + s ] c ArcSinh[1 - s] + c ArcSinh[s]) We now determine the least square function of the parameters a and b in our function x(t) = a + b t + c t^2. We seek to minimize ss(a, b,c) . ss[a_,b_,c_] = Sum[(zz[data[[i,1]]] - data[[i,2]])^2, {i,1,Length[data]}]//N; We use Mathematica's FindMinimum command to determine the values of a and b which fit the z(s) data best. sol = FindMinimum[ss[a,b,c],{a,3},{b,2},{c,4}] -25 {1.82481 10 , {a -> 0.839822, b -> 1.73254, c -> -1.93651}} We create our best fitting model function. xb[t_] = x[t]/.sol[[2]] 2 0.839822 + 1.73254 t - 1.93651 t So how good is our model for x(t), i.e. xb(t)? We put it into our z(s) function and plot the theoretical z's against the observed ones. z[s_] = G Integrate[((s-t)^2 + 1)^(-3/2) xb[t],{t,0,1}] -11 6.67 10 ( 2 3 1.04378 + 2.82923 s - 3.66905 s + 1.93651 s --------------------------------------------- - 2 Sqrt[2 - 2 s + s ] 2 3 -1.73254 + 1.09668 s - 1.73254 s + 1.93651 s ---------------------------------------------- - 2 Sqrt[1 + s ] 1.93651 ArcSinh[1 - s] - 1.93651 ArcSinh[s]) zplot = Plot[z[s],{s,0,1}] -Graphics- Not all that bad!!!!! Show[ysp,zplot] -Graphics- ISSUES IN SOLUTION Students do have difficulty with finding and adding up the vertical components of the forces. There are a number of different concepts running through this problem and students need to keep straight where they are in the analysis. Summarizing what they have accomplished at each stage is important before they proceed.