An object's center of mass provides the point where the object can balance perfectly. In physics and engineering, the center of mass gives great insight into how an object will behave in certain instances, and gives a point of reference for where forces on an object will act.
Mathematically, we can find the center of mass of a plate of uniform density, by finding the area enclosed by a set of curves.
I use two equations to define my heart shape, and define them over the interval \([-1, 1]\)
\[f(x) = |x|^\frac{2}{3}\ + \sqrt[]{1 - x ^2}\]
\[g(x) = |x|^\frac{2}{3}\ - \sqrt[]{1 - x ^2}\]
Below I give the plot for my heart, and the rendering of its STL file.
To identify the moment of this object with respect to the y-axis over the x-interval \([a, b]\), we use the following formula. \(\rho\) represents the uniform density of the 'material' that consists of our plate.
\[M_y = \rho \int_{a}^{b} x f(x) dx\]
We do the same thing for the moment with respect to the x-axis, with the y-interval \([c, d]\), where \(\rho\) once again represents the uniform density of our plate 'material'.
\[M_x = \rho \int_{c}^{d} y f(y) dy\]
Both of the above equations assume that we do not need to separate intervals of x or y. However we can, because moments are additive! We will actually have to do this for our functions above for the heart, since it does not have continuity due to the absolute value, and because it comprises of two separate equations.
To find the mass of our plate, we simply integrate over the region defined by the curve \(f(x)\), over the interval \([a,b]\). And \(\rho\) still represents uniform density.
\[Mass = \rho \int_{a}^{b} f(x) dx\]
Though in my case, I have two curves, so it yields the combined expression below.
\[Mass = \rho \int_{a}^{b} (f(x) + g(x) )dx\]
Thankfully, the sum of my functions yields the wonderfully simple expression below. Since I can only integrate the interval \([0,1]\) I have to multiply this expression by two. Likewise, the sum of the \(x^\frac{2}{3}\) term causes me to multiply by two again, and our radical terms sum to zero!
\[Mass = 4 \rho \int_{0}^{1} x^\frac{2}{3} dx \]
\[Mass = 4 \rho (\frac{3}{5} ( x^\frac{5}{3} \Big|_{0}^{1}))\]
\[Mass = \frac{12}{5} \rho\]
We can combine these expressions together to obtain the center of mass. This boils down to the moment with respect to the x-axis divided by the mass, and the moment with respect to the y-axis divided by the mass.
\[C.O.M. = (\frac{M_x}{Mass}, \frac{M_y}{Mass} )\]
We will take advantage of the additive property of moments and begin with finding the moment with respect to the y-axis. Since our first curve \(f(x)\) has an absolute value, we can take advantage of the symmetry of this function, and integrate only over the interval \([0,1]\) ! So, the value we obtain on the interval \([-1, 0]\) has the negative value of what we obtain on \([0,1]\). This fact will also hold for our second function \(g(x)\).
\[M_y{_1}{_{rightside}} = \rho \int_{0}^{1} x (x^\frac{2}{3} +\sqrt[]{1 - x^2}) dx\]
\[M_y{_1}{_{rightside}} = \rho (\frac{3}{8} x^\frac{8}{3} + \frac{1}{3} x ^\frac{3}{2} \Big|_{0}^{1})\]
\[M_y{_1}{_{rightside}} = \rho \frac{17}{24} \]
Since we have symmetry about the y-axis...
\[M_y{_1}{_{leftside}} = -M_y{_1}{_{rightside}} \]
So we sum the left side with the right side and we get for the first equation \(f(x)\)
\[M_y{_1} = 0\]
We have to repeat this process for \(g(x)\)
\[M_y{_2}{_{rightside}} = \rho \int_{0}^{1} x (x^\frac{2}{3} -\sqrt[]{1 - x^2}) dx\]
\[M_y{_2}{_{rightside}} = \rho (\frac{3}{8} x^\frac{8}{3} - \frac{1}{3} x ^\frac{3}{2} \Big|_{0}^{1})\]
\[M_y{_2}{_{rightside}} = \rho \frac{1}{24} \]
Once again we take advantage of the symmetry about the y-axis for this function, and we see that the left side has the negative value of the right side.
\[M_y{_2}{_{leftside}} = -M_y{_2}{_{rightside}} \]
\[M_y{_2} = 0\]
To obtain our total moment about the y-axis we have to sum together \(M_y{_1}\) and \(M_y{_2}\).
\[M_y = M_y{_1} + M_y{_2} \\ M_y = 0 + 0 \\ M_y = 0\]
Now to solve for the moment about the x-axis, we have to solve our functions of x in terms of y... This yields the awful expressions below. Though, as seen in the curve above, \(f(x)\) is not one-to-one in terms of y, and \(y(0) = y(1) = 1\). While this should have an analytical solution, I refuse, refuse to compute it. Here is Wolfram Alpha's attempt.
Below I give the code I use to numerically approximate the moment with respect to the x-axis for \(f(x)\) and \(g(x)\). Below I give the plot for y in terms of x.
My scheme to find the solution is as follows:
- Separate curve about x = 0.
- Find the array that contains y = 1 to y = max(y)
- Obtain two dictionaries which map x-values to the range of y=1 to y = max(y)
- Find the element-wise difference between the two sets of keys from the two dictionaries
- Multiply the array of differences by the range of y-values to give the moment
- Multiply the product of the y-range by an appropriately small dy value (my domain has 10,000 steps between 0 and 1, so \(\displaystyle dy = \frac{1}{2500}\)
- Take the sum of this array and multiply by two to give the other half of the area.
- Repeat for \(g(x)\)
- Add together the sum of \(g(x)\) in terms of y and \(f(x)\) in terms of y to give the total moment with respect to the x-axis
\[X Axis Moment_f = 1.349 \\ X Axis Moment_g = 0.0835 \\ M_x = 1.4335\]
The center of mass is as follows:
\[COM = (\frac{0}{\frac{12}{5}} , \frac{1.4335}{\frac{12}{5}})\]
\[COM = (0, 0.5973)\]
My code is given below:
import numpy as np import matplotlib.pyplot as plt def give_single_y_for_f(x): return (x**(2/3) + np.sqrt(1-x**2)) def give_single_y_for_g(x): return (x**(2/3) - np.sqrt(1-x**2)) def find_x_moment(x_values, y_values, which ="f"): positive_domain = x_values[x_values>0] ymin = y_values.min() ymax = y_values.max() upper_dict = {} lower_dict = {} y_prev = -999 for x in positive_domain: if which == 'f': y=give_single_y_for_f(x) else: y = give_single_y_for_g(x) if y>y_prev: lower_dict[x] = y y_prev = y else: upper_dict[x] = y y_prev = y upper_keys_as_array = np.array(list(upper_dict.keys())) lower_keys_as_array = np.array(list(lower_dict.keys())) if len(upper_keys_as_array)!= len(lower_keys_as_array) and which =="f": _min = upper_keys_as_array.min() _max = upper_keys_as_array.max() upper_keys_as_array = np.linspace(_min, _max, len(lower_keys_as_array)) difference_array = upper_keys_as_array-lower_keys_as_array else: difference_array = lower_keys_as_array y_range = np.array(list(lower_dict.values())) moment_array = y_range*difference_array dy = 1/len(y_range) y_dy = dy*moment_array total_moment = sum(y_dy) if which =="f": return 2*total_moment else: return 2*total_moment pos_domain = np.linspace(-1,1, 10000) pos_two_thirds = np.sqrt((np.abs(pos_domain))**3) reflect_two_thirds = -pos_two_thirds pos_one_minus = np.sqrt(1-(np.abs(pos_domain)**2)) pos_half = pos_two_thirds+pos_one_minus other_pos_half = pos_two_thirds - pos_one_minus print(find_x_moment(pos_domain, pos_half)) print(find_x_moment(pos_domain, other_pos_half, which = "g")) moment_one = find_x_moment(pos_domain, other_pos_half, which = "g") moment_two = find_x_moment(pos_domain, pos_half) the_moment_sum = moment_one+moment_two print(F"\n\tThe total moment about the x-axis is {round(the_moment_sum,5)}") plt.plot(pos_domain, pos_half) plt.plot(pos_domain, other_pos_half) plt.legend(["f(x)", "g(x)"]) plt.title("My Lame Heart") plt.ylim([-1.5, 1.5]) plt.show() plt.plot(pos_half, pos_domain) plt.title("f(y) = pain") plt.show() plt.plot(other_pos_half, pos_domain) plt.title("g(y) = more pain") plt.show()
Comments
Post a Comment