Original file (1,800 × 1,200 pixels, file size: 19.36 MB, MIME type: image/gif, looped, 201 frames, 16 s)
Note: Due to technical limitations, thumbnails of high resolution GIF images such as this one will not be animated.
This is a file from the
Wikimedia Commons. Information from its
description page there is shown below. Commons is a freely licensed media file repository. You can help. |
DescriptionPhase portrait of damped oscillator, with increasing damping strength.gif |
English: ```python
import numpy as np import matplotlib.pyplot as plt from math import isclose from numpy import linalg as LA import matplotlib.cm as cm def plot_circle(v_1, v_2, ax, **kwargs): angles = np.linspace(0, 2*np.pi, 100) points = v_1[:,np.newaxis] * np.cos(angles) + v_2[:,np.newaxis] * np.sin(angles) ax.plot(points[0,:], points[1,:], **kwargs) def plot_vector_field(A, xmin=-5, xmax=5, ymin=-5, ymax=5, title=""): axx, axy = A[0,0], A[0,1] ayx, ayy = A[1,0], A[1,1] det = axx * ayy - axy * ayx tr = axx + ayy eigen_vals, eigen_vects = LA.eig(A) is_critical = abs(eigen_vals[0] - eigen_vals[1]) / abs(eigen_vals[0]) < 1e-2 delta = tr**2 - 4*det is_rotational = delta <= 0 and not is_critical # Initialize plotting object fig, axes = plt.subplot_mosaic("133;233", figsize=(18,12)) colormap=cm.viridis # pole-zero plot ax = axes['1'] ax.scatter(eigen_vals[0].real, eigen_vals[0].imag, color=colormap(eigen_vals[0].real)) ax.scatter(eigen_vals[1].real, eigen_vals[1].imag, color=colormap(eigen_vals[1].real)) r = np.sqrt(abs(eigen_vals[0] * eigen_vals[1])) plot_circle(np.array([r, 0]), np.array([0, r]), ax, color='k', alpha=0.3) ax.plot([xmin, xmax], np.zeros(2), color='k', alpha=0.2) ax.plot(np.zeros(2), [ymin, ymax], color='k', alpha=0.2) ax.set_aspect('equal') ax.set_xlim([-2, 2]) ax.set_ylim([-2, 2]) ax.set_xlabel('Real') ax.set_ylabel('Imag') ax.set_title('pole-zero plot') # stability plot ax = axes['2'] xs = np.linspace(xmin, xmax, 100) ys = xs**2 / 4 ax.plot(xs, ys) ax.scatter(tr, det, color='red') ax.plot([xmin, xmax], np.zeros(2), color='k', alpha=0.2) ax.plot(np.zeros(2), [ymin, ymax], color='k', alpha=0.2) ax.set_aspect('equal') ax.set_xlim([-4,2]) ax.set_ylim([-1, 5]) ax.set_xlabel('Tr(A)') ax.set_ylabel('Det(A)') ax.set_title('stability plot') # vector field plot ax = axes['3'] x, y = np.meshgrid(np.linspace(xmin, xmax, 10), np.linspace(ymin, ymax, 10)) vx = axx * x + axy * y vy = ayx * x + ayy * y ax.quiver(x,y, vx, vy, units='xy', scale=6, color='g', headwidth=3, width=0.04) # Plot the circle, or fast and slow manifolds if is_rotational: v_1 = np.array(eigen_vects[:,0].real) v_2 = np.array(eigen_vects[:,0].imag) # normalize radius = (xmax - xmin) / 4 norm = max(np.linalg.norm(v_1), np.linalg.norm(v_2)) / radius v_1 /= norm v_2 /= norm plot_circle(v_1, v_2, ax, color=colormap(eigen_vals[0].real)) elif is_critical: v_1 = eigen_vects[:,0] length = (xmax - xmin) * 2 lengths = np.linspace(-length, length, 100) points = v_1[:,np.newaxis] * lengths ax.plot(points[0,:], points[1,:], color=colormap(eigen_vals[0])) else: v_1 = eigen_vects[:,0] v_2 = eigen_vects[:,1] length = (xmax - xmin) * 2 lengths = np.linspace(-length, length, 100) points = v_1[:,np.newaxis] * lengths ax.plot(points[0,:], points[1,:], color=colormap(eigen_vals[0])) points = v_2[:,np.newaxis] * lengths ax.plot(points[0,:], points[1,:], color=colormap(eigen_vals[1])) ax.plot([xmin, xmax], np.zeros(2), color='k', alpha=0.2) ax.plot(np.zeros(2), [ymin, ymax], color='k', alpha=0.2) ax.set_aspect('equal') ax.set_xlim([xmin, xmax]) ax.set_ylim([ymin, ymax]) ax.set_xlabel('$x$') ax.set_ylabel('$\\dot{x}$') ax.set_title('phase portrait') fig.suptitle(title) fig.tight_layout() return fig import tempfile import os import imageio plt.rc('figure', titlesize=16) with tempfile.TemporaryDirectory() as temp_dir: n_frames = 201 omegas = [1.0] * n_frames gammas = (1-np.cos(np.linspace(0, np.pi, n_frames//2)))/2 gammas = list(gammas) + [gammas[-1]] + list(gammas + 1) for i in range(n_frames): omega = omegas[i] gamma = gammas[i] operator_A = np.array([[0, 1], [-omega**2, -2*gamma]]) fig = plot_vector_field(operator_A, title=f"$\\ddot x + 2\\gamma \\dot x + \\omega^2x = 0$,\n$\\omega = {omega:1.1f}$, $\\gamma = {gamma:0.3f}$") filename = os.path.join(temp_dir, f"plot_{i:03d}.png") fig.savefig(filename) plt.close(fig) # Compile images into GIF fps = 12 images = [] for i in range(n_frames): filename = os.path.join(temp_dir, f"plot_{i:03d}.png") images.append(imageio.v2.imread(filename)) imageio.mimsave(f"phase_portrait_omega_{omega:1.1f}.gif", images, duration=1/fps)``` |
Date | |
Source | Own work |
Author | Cosmia Nebula |
Click on a date/time to view the file as it appeared at that time.
Date/Time | Thumbnail | Dimensions | User | Comment | |
---|---|---|---|---|---|
current | 06:24, 5 April 2023 | 1,800 × 1,200 (19.36 MB) | Cosmia Nebula | Uploaded own work with UploadWizard |
The following other wikis use this file: