Skip to main content

Center of Mass of a 'Heart'

 

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
 Using this method, I obtained:
\[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()















Labels:

Comments

Post a Comment

Popular posts from this blog

Do Over: Integration Over a Region in a Plane

Throughout the semester we have covered a variety of topics and how their mathematical orientation applies to real world scenarios. One topic we discussed, and I would like to revisit, is integration over a region in a plane which involves calculating a double integral. Integrating functions of two variables allows us to calculate the volume under the function in a 3D space. You can see a more in depth description and my previous example in my blog post, https://ukyma391.blogspot.com/2021/09/integration-for-over-regions-in-plane_27.html . I want to revisit this topic because in my previous attempt my volume calculations were incorrect, and my print lacked structural stability. I believed this print and calculation was the topic I could most improve on and wanted to give it another chance. What needed Improvement? The function used previously was f(x) = cos(xy) bounded on [-3,3] x [-1,3]. After solving for the estimated and actual volume, it was difficult to represent in a print...
Post a Comment
Read more

Minimal Surfaces

Minimum surfaces can be described in many equivalent ways. Today, we are going to focus on minimum surfaces by defining it using curvature. A surface is a minimum surface if and only if the mean curvature at every point is zero. This means that every point on the surface is a saddle point with equal and opposite curvature allowing the smallest surface area possible to form. Curvature helps define a minimal surface by looking at the normal vector. For a surface in R 3 , there is a tangent plane at each point. At each point in the surface, there is a normal vector perpendicular to the tangent plane. Then, we can intersect any plane that contains the normal vector with the surface to get a curve. Therefore, the mean curvature of a surface is defined by the following equation. Where theta is an angle from a starting plane that contains the normal vector. For this week’s project, we will be demonstrating minimum surfaces with a frame and soap bubbles! How It Works Minimum surfac...
Post a Comment
Read more

Ruled Surfaces : Trefoil

Ruled Surfaces : Trefoil A ruled surface is a surface that consists straight lines, called rulings, which lie upon the surface. These surfaces are formed of a set of points that are "swept" by a straight line. This is relatively intuitive once you see a good visual, but can be a bit abstract without that concrete example. A very basic example of a ruled surface is a cylinder; if we have a straight line and move it in a circle we create a cylinder made entirely of straight line. Note that the surface will only be a cylinder if all the lines are parallel. If the lines are not parallel we can create hyperboloids and cones depending on how much we have rotated. The rotation we are describing here is not a simple turning action, but more of a twisting motion—less like rotating a can by turning it and more like wringing out a washcloth by twisting it. Specifically, a cylinder is essentially two circles connected by rulings, if we keep one of the circles...
Post a Comment
Read more