Happy Valentine’s Linkage#
In this notebook, you step through the process of defining the kinematics of a four-bar linkage that can draw a heart. Kinematics is the study of the geometry of motion. In this notebook, you predict and draw the geometry of three moving links. To create this solution you will:
solve a series of nonlinear equations using
fsolveuse solutions to create 2D arrays that vary in time and location
plot and animate the motion of the four-bar linkage
When you accomplish these steps you will learn:
How to set up and solve nonlinear equations
How to describe position and orientation
How to use NumPy, Scipy, and Matplotlib to create HTML animations
Why four-bar linkages are so cool
Background#
The four-bar linkages consist of 3 moving parts connected to a stationary support. Depending upon the desired output motion, you can vary the lengths of each moving arm, \(l_1,~l_2,~and,l_3\), and the support locations, \(d_x~and~d_y\). A four-bar linkage is an amazing 1-degree-of-freedom system. If one angle is fixed, all of the positions and angles have to be fixed too. Using some trigonometry, you can arrive at these two constraint equations that relate \(\theta_1,~\theta_2,~and~\theta_3\):
\(l_1\sin\theta_1+l_2\sin\theta_2-l_3\sin\theta_3 -d_y = 0\)
\(l_1\cos\theta_1+l_2\cos\theta_2-l_3\cos\theta_3 -d_x = 0\)
If you have one of the angles, e.g. \(\theta_1\), you use equations 1 and 2
to solve for the other two angles, \(\theta_2~and\theta_3\). Here you can
create a function and
use fsolve. The function input is a vector with two values and the output is a
vector with two values.
\(\bar{f}(\bar{x})= \left[\begin{array}{c} f_1(\theta_2,~\theta_3) \\ f_2(\theta_2,~\theta_3)\end{array}\right]=\left[\begin{array}{c} l_1\sin\theta_1+l_2\sin\theta_2-l_3\sin\theta_3 -d_y\\ l_1\cos\theta_1+l_2\cos\theta_2-l_3\cos\theta_3 -d_x \end{array}\right]\)
Defining your system#
The heart-drawing linkage system has the following properties:
link 1: \(l_1 = 1~m\)
link 2: \(l_2 = 1~m\)
link 3: \(l_3 = 1~m\)
support: \(d_x=0.95~m~and~d_y=0~m\)
The constraint function is defined below as Fbar, a function of
\(\theta_1\) and an array of \([\theta_2,~\theta_3]\) as such,
l1 = 1
l2 = 1
l3 = 1
a1 = np.pi/2
dy = 0
dx = 0.95
Fbar = lambda a1,x: np.array([l1*np.sin(a1)+l2*np.sin(x[0])-l3*np.sin(x[1])-dy,
l1*np.cos(a1)+l2*np.cos(x[0])-l3*np.cos(x[1])-dx])
Solve for one configuration#
Now, solve Fbar using
scipy.optimize.fsolve.
The inputs are a function, Fbar, and an initial guess, x0. You have
to use lambda again to set the angle \(\theta_1\) as a1 = np.pi/2
a1 = np.pi/2
x0 = np.array([0,np.pi/2])
xsol = fsolve(lambda x: Fbar(a1, x), x0)
The configuration for \(\theta_1=\frac{\pi}{2}\) is now saved in xsol.
To look at the system in this state, define the positions of each hinge
and the center of link 2 as such,
x- and y-locations of hinges:
\(rx = \left[\begin{array}~0\\l_1\cos(\theta_1)\\l_1\cos(\theta_1)+l_2\cos(\theta_2)\\ l_1\cos(\theta_1) + l_2\cos(\theta_2)-l_3\cos(\theta_3)\end{array}\right]\)
\(ry = \left[\begin{array}~0\\l_1\sin(\theta_1)\\l_1\sin(\theta_1)+l_2\sin(\theta_2)\\ l_1\sin(\theta_1)+l_2\sin(\theta_2)-l_3\sin(\theta_3)\end{array}\right]\)
x- and y-location of point P:
\(rx = \left[\begin{array}~l_1\cos(\theta_1)+l_2\cos(\theta_2)\end{array}\right]\)
\(ry = \left[\begin{array}~l_1\sin(\theta_1)+l_2\sin(\theta_2)\end{array}\right]\)
In the Python cell, you define rx and ry to draw the three moving
links and rpx and rpy that define point P as such,
rx = np.array([0,
l1*np.cos(a1),
l1*np.cos(a1)+l2*np.cos(xsol[0]),
l1*np.cos(a1)+l2*np.cos(xsol[0])-l3*np.cos(xsol[1])])
ry = np.array([0,
l1*np.sin(a1),
l1*np.sin(a1)+l2*np.sin(xsol[0]),
l1*np.sin(a1)+l2*np.sin(xsol[0])-l3*np.sin(xsol[1])])
rpx = l1*np.cos(a1)+l2/2*np.cos(xsol[0])
rpy = l1*np.sin(a1)+l2/2*np.sin(xsol[0])
plt.plot(rx,ry,'o-')
plt.plot(rpx,rpy,'rs')
#plt.axis([-0.5, 1.5, -0.6, 0.6])
[<matplotlib.lines.Line2D at 0x7f36377a1e90>]
Solve for the whole cycle#
You have verified the solution works for one angle, a1. Now, you can
create an array of a1 and solve for the angles at each configuration
as you rotate link 1.
define
a1as alinspaceinitialize
a2anda3as zeros witha1.shapeuse a for-loop to
fsolveeach anglea2anda3givena1save the values of
a2anda3in each step
a1 = np.linspace(-np.pi/2,3*np.pi/2,500)
a2 = np.zeros(a1.shape) # initialize
a3 = np.zeros(a1.shape) # initialize
for i, a in enumerate(a1):
xsol = fsolve(lambda x: Fbar(a,x), xsol) # solve
a2[i] = xsol[0] # save value for a2
a3[i] = xsol[1] # save value for a3
Verify cycle solution#
Now, you can see how the functions, \(\theta_2=f_{\theta_2}(\theta_1)\) and \(\theta_3=f_{\theta_3}(\theta_1)\). This is another verification step, to see if there are any discontinuities in the solution that could cause problems.
plt.plot(a1*180/np.pi, a2*180/np.pi, label = r'$\theta_2$')
plt.plot(a1*180/np.pi, a3*180/np.pi, label = r'$\theta_3$')
plt.legend()
plt.xlabel(r'input: $\theta_1$ (degrees)')
plt.ylabel('output: angles (degrees)')
Text(0, 0.5, 'output: angles (degrees)')
Time to animate#
Now, you are ready to animate. You have the angles for the entire cycle
of motion for the four-bar linkage. Now, you need to import
matplotlib.animation
and
IPython.display.HTML.
Import functions#
from matplotlib import animation
from IPython.display import HTML
These two functions allow you:
create an animation
display it in a browser
Define lines and paths#
Then, use the solutions for a1, a2, and a3 to plot the hinge
locations and point P as such
rx = np.array([np.zeros(a1.shape),
l1*np.cos(a1),
l1*np.cos(a1)+l2*np.cos(a2),
l1*np.cos(a1)+l2*np.cos(a2)-l3*np.cos(a3)])
ry = np.array([np.zeros(a1.shape),
l1*np.sin(a1),
l1*np.sin(a1)+l2*np.sin(a2),
l1*np.sin(a1)+l2*np.sin(a2)-l3*np.sin(a3)])
rpx = l1*np.cos(a1)+l2/2*np.cos(a2)
rpy = l1*np.sin(a1)+l2/2*np.sin(a2)
plt.plot(rpx,rpy)
[<matplotlib.lines.Line2D at 0x7f3637184510>]
Set up figure and axes#
Plotting just the solution for point P’s path, you see the heart shape is upside-down. In the animation, you can reverse the y-axis to flip the drawing. Here, you set up the figure to create the animation:
ax: axis for plotting the linesline1: lines that draw the three moving linksrxandryline2: line that updates the heart drawing as links move along paths
fig, ax = plt.subplots()
ax.set_ylim((1.5, -1.5))
ax.set_xlabel('x-position (m)')
ax.set_ylabel('y-position (m)')
ax.set_aspect('equal')
line1, = ax.plot([], [],'bo-')
line2, = ax.plot([], [],'r')
ax.plot(rx[1,:],ry[1,:],'g--', alpha=0.5)
ax.plot(rx[2,:],ry[2,:],'g--', alpha=0.5)
[<matplotlib.lines.Line2D at 0x7f3636bef250>]
Define initializing and animation functions#
Create an initializing (init) function that clears the previous
lines
def init():
line1.set_data([], [])
line2.set_data([], [])
return (line1, line2, )
Create an animating (animate) function that updates the lines. The
links should be drawn for a given column, i.e. line1 is defined as
rx[:, i] and ry[:, i], but you want the path of the heart up to the
current frame i.e. line2 is defined as rpx[:i] and rpy[:i].
def animate(i):
'''function that updates the line and marker data
arguments:
----------
i: index of timestep
outputs:
--------
line: the line object plotted in the above ax.plot(...)
'''
line1.set_data(rx[:, i], ry[:,i])
line2.set_data(rpx[:i], rpy[:i])
return (line1, line2 )
Create the animation and display#
Create an animation (anim) variable using the animation.FuncAnimation
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=range(0,len(a1)), interval=20,
blit=True)
HTML(anim.to_html5_video())
MovieWriter stderr:
[out#0/ipod @ 0x558585dae680] Error writing trailer: Immediate exit requested
[out#0/ipod @ 0x558585dae680] Error closing file: Immediate exit requested
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:224, in AbstractMovieWriter.saving(self, fig, outfile, dpi, *args, **kwargs)
223 try:
--> 224 yield self
225 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1126, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1125 frame_number += 1
-> 1126 writer.grab_frame(**savefig_kwargs)
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:352, in MovieWriter.grab_frame(self, **savefig_kwargs)
351 # Save the figure data to the sink, using the frame format and dpi.
--> 352 self.fig.savefig(self._proc.stdin, format=self.frame_format,
353 dpi=self.dpi, **savefig_kwargs)
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/figure.py:3490, in Figure.savefig(self, fname, transparent, **kwargs)
3489 _recursively_make_axes_transparent(stack, ax)
-> 3490 self.canvas.print_figure(fname, **kwargs)
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/backend_bases.py:2186, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2185 with cbook._setattr_cm(self.figure, dpi=dpi):
-> 2186 result = print_method(
2187 filename,
2188 facecolor=facecolor,
2189 edgecolor=edgecolor,
2190 orientation=orientation,
2191 bbox_inches_restore=_bbox_inches_restore,
2192 **kwargs)
2193 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/backend_bases.py:2042, in FigureCanvasBase._switch_canvas_and_return_print_method.<locals>.<lambda>(*args, **kwargs)
2041 skip = optional_kws - {*inspect.signature(meth).parameters}
-> 2042 print_method = functools.wraps(meth)(lambda *args, **kwargs: meth(
2043 *args, **{k: v for k, v in kwargs.items() if k not in skip}))
2044 else: # Let third-parties do as they see fit.
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:417, in FigureCanvasAgg.print_raw(self, filename_or_obj, metadata)
416 raise ValueError("metadata not supported for raw/rgba")
--> 417 FigureCanvasAgg.draw(self)
418 renderer = self.get_renderer()
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:382, in FigureCanvasAgg.draw(self)
380 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
381 else nullcontext()):
--> 382 self.figure.draw(self.renderer)
383 # A GUI class may be need to update a window using this draw, so
384 # don't forget to call the superclass.
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/artist.py:94, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
92 @wraps(draw)
93 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 94 result = draw(artist, renderer, *args, **kwargs)
95 if renderer._rasterizing:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/figure.py:3257, in Figure.draw(self, renderer)
3256 self.patch.draw(renderer)
-> 3257 mimage._draw_list_compositing_images(
3258 renderer, self, artists, self.suppressComposite)
3260 renderer.close_group('figure')
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/axes/_base.py:3226, in _AxesBase.draw(self, renderer)
3224 _draw_rasterized(self.get_figure(root=True), artists_rasterized, renderer)
-> 3226 mimage._draw_list_compositing_images(
3227 renderer, self, artists, self.get_figure(root=True).suppressComposite)
3229 renderer.close_group('axes')
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/image.py:134, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
133 for a in artists:
--> 134 a.draw(renderer)
135 else:
136 # Composite any adjacent images together
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/axis.py:1408, in Axis.draw(self, renderer)
1407 for tick in ticks_to_draw:
-> 1408 tick.draw(renderer)
1410 # Shift label away from axes to avoid overlapping ticklabels.
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/axis.py:276, in Tick.draw(self, renderer)
274 for artist in [self.gridline, self.tick1line, self.tick2line,
275 self.label1, self.label2]:
--> 276 artist.draw(renderer)
277 renderer.close_group(self.__name__)
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/artist.py:71, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 renderer.start_filter()
---> 71 return draw(artist, renderer)
72 finally:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/lines.py:821, in Line2D.draw(self, renderer)
820 gc.set_dashes(*self._dash_pattern)
--> 821 renderer.draw_path(gc, tpath, affine.frozen())
822 gc.restore()
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/backends/backend_agg.py:130, in RendererAgg.draw_path(self, gc, path, transform, rgbFace)
129 try:
--> 130 self._renderer.draw_path(gc, path, transform, rgbFace)
131 except OverflowError:
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/backend_bases.py:800, in GraphicsContextBase.get_snap(self)
798 return self._gid
--> 800 def get_snap(self):
801 """
802 Return the snap setting, which can be:
803
(...) 807 round to the nearest pixel center
808 """
KeyboardInterrupt:
During handling of the above exception, another exception occurred:
CalledProcessError Traceback (most recent call last)
Cell In[14], line 1
----> 1 HTML(anim.to_html5_video())
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1306, in Animation.to_html5_video(self, embed_limit)
1302 Writer = writers[mpl.rcParams['animation.writer']]
1303 writer = Writer(codec='h264',
1304 bitrate=mpl.rcParams['animation.bitrate'],
1305 fps=1000. / self._interval)
-> 1306 self.save(str(path), writer=writer)
1307 # Now open and base64 encode.
1308 vid64 = base64.encodebytes(path.read_bytes())
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:1098, in Animation.save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1093 return a * np.array([r, g, b]) + 1 - a
1095 # canvas._is_saving = True makes the draw_event animation-starting
1096 # callback a no-op; canvas.manager = None prevents resizing the GUI
1097 # widget (both are likewise done in savefig()).
-> 1098 with (writer.saving(self._fig, filename, dpi),
1099 cbook._setattr_cm(self._fig.canvas, _is_saving=True, manager=None)):
1100 if not writer._supports_transparency():
1101 facecolor = savefig_kwargs.get('facecolor',
1102 mpl.rcParams['savefig.facecolor'])
File ~/miniconda3/envs/jupbook01/lib/python3.11/contextlib.py:158, in _GeneratorContextManager.__exit__(self, typ, value, traceback)
156 value = typ()
157 try:
--> 158 self.gen.throw(typ, value, traceback)
159 except StopIteration as exc:
160 # Suppress StopIteration *unless* it's the same exception that
161 # was passed to throw(). This prevents a StopIteration
162 # raised inside the "with" statement from being suppressed.
163 return exc is not value
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:226, in AbstractMovieWriter.saving(self, fig, outfile, dpi, *args, **kwargs)
224 yield self
225 finally:
--> 226 self.finish()
File ~/miniconda3/envs/jupbook01/lib/python3.11/site-packages/matplotlib/animation.py:341, in MovieWriter.finish(self)
337 _log.log(
338 logging.WARNING if self._proc.returncode else logging.DEBUG,
339 "MovieWriter stderr:\n%s", err)
340 if self._proc.returncode:
--> 341 raise subprocess.CalledProcessError(
342 self._proc.returncode, self._proc.args, out, err)
CalledProcessError: Command '['ffmpeg', '-f', 'rawvideo', '-vcodec', 'rawvideo', '-s', '640x480', '-pix_fmt', 'rgba', '-framerate', '50.0', '-loglevel', 'error', '-i', 'pipe:', '-vcodec', 'h264', '-pix_fmt', 'yuv420p', '-y', '/tmp/tmp49xc_fie/temp.m4v']' returned non-zero exit status 255.
Wrapping up#
There you have it, a heart-drawing four-bar linkage. This solution required arrays, nonlinear solutions, and some trigonometry to build locations of hinges and links.
Try changing the lengths the different links or changing the fixed positions. What else can you draw?