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

Do Over: Ruled Surfaces

Why to choose this project to repeat For the do over project, I would like to choose the ruled surfaces. I don't think my last project was creative, and the 3D printed effect was not very satisfactory. In the previous attempts, all the lines are connected between a straight line and a circle. This connection structure is relatively uncomplicated. The printed model has too many lines, resulting in too dense line arrangement. The gaps between lines are too small, and the final effect is that all the lines are connected into a curved surface, which is far from the effect I expected. What to be improved In this do over project, I would like to improve in two aspects. Firstly, a different ruled surface is chosen. In the previous model, one curve is a unit circle on the \(x-y\) plane, and the ruled surface is a right circular conoid. In this do over project, it is replaced by two border lines. Each borderline is in the shape of an isosceles right triangl...
Post a Comment
Read more