Python科学计算中arctan函数的实现方法与应用
作者:麦兜 更新时间:2023-06-28 21:00:01 阅读 934
散点图
import numpy as np
import matplotlib.pyplot as plt
countries = ['Brazil', 'Madagascar', 'S. Korea', 'United States',
'Ethiopia', 'Pakistan', 'China', 'Belize']
# Birth rate per 1000 population
birth_rate = [16.4, 33.5, 9.5, 14.2, 38.6, 30.2, 13.5, 23.0]
# Life expectancy at birth, years
life_expectancy = [73.7, 64.3, 81.3, 78.8, 63.0, 66.4, 75.2, 73.7]
# Per person income fixed to US Dollars in 2000
GDP = np.array([4800, 240, 16700, 37700, 230, 670, 2640, 3490])
fig = plt.figure()
ax = fig.add_subplot(111)
# Some random colours:
colours = range(len(countries))
ax.scatter(birth_rate, life_expectancy, c=colours, s=GDP/20)
ax.set_xlim(5, 45)
ax.set_ylim(60, 85)
ax.set_xlabel('Birth rate per 1000 population')
ax.set_ylabel('Life expectancy at birth (years)')
plt.show()
要给气泡标上标签，在显示图形之前添加以下代码:
offset = 1
for x, y, s, country in zip(birth_rate, life_expectancy, GDP, countries):
ax.text(x+offset, y, country, va='center')
曲线图
数据eg7-marriage-ages.txt
Median age at First Marriage, 1890-2010
Source: U.S. Bureau of the Census, www.census.gov
Year Men Women
1890 26.1 22.0
1900 25.9 21.9
1910 25.1 21.6
1920 24.6 21.2
1930 24.3 21.3
1940 24.3 21.5
1950 22.8 20.3
1960 22.8 20.3
1970 23.2 20.8
1980 24.7 22.0
1990 26.1 23.9
2000 26.8 25.1
2010 28.2 26.1import numpy as np
import matplotlib.pyplot as plt
year, age_m, age_f = np.loadtxt('eg7-marriage-ages.txt', unpack=True,skiprows=3)
fig = plt.figure()
ax = fig.add_subplot(111)
# Plot ages with male or female symbols as markers
ax.plot(year, age_m, marker='$\u2642#39;, markersize=14, c='blue', lw=2,
mfc='blue', mec='blue')
ax.plot(year, age_f, marker='$\u2640#39;, markersize=14, c='magenta', lw=2,
mfc='magenta', mec='magenta')
ax.grid()
ax.set_xlabel('Year')
ax.set_ylabel('Age')
ax.set_title('Median age at first marriage in the US, 1890 - 2010')
plt.show()
以下程序将这些数据绘制在一组轴上，每个轴具有不同的线型。
import matplotlib.pyplot as plt
import numpy as np
fig = plt.figure()
ax = fig.add_subplot(111)
cities = ['Boston', 'Houston', 'Detroit', 'San Jose', 'Phoenix']
# line styles: solid, dashes, dots, dash-dots, and dot-dot-dash
linestyles = [{'ls': '-'}, {'ls': '--'}, {'ls': ':'}, {'ls': '-.'},
{'dashes': [2, 4, 2, 4, 8, 4]}]
for i, city in enumerate(cities):
filename = '{}.tsv'.format(city.lower()).replace(' ', '_')
yr, pop = np.loadtxt(filename, unpack=True)
line, = ax.plot(yr, pop/1.e6, label=city, color='k', **linestyles[i])
ax.legend(loc='upper left')
ax.set_xlim(1800, 2020)
ax.set_xlabel('Year')
ax.set_ylabel('Population (millions)')
plt.show()
请注意如何使用城市名称来推断相应的文件名。
指数衰减的寿命
下面的程序绘制了用寿命标记的描述的指数衰减，(nτ for n=0,1，⋯)这样，在每个寿命之后，y的值下降了一个因子e。
import numpy as np
import matplotlib.pyplot as plt
# Initial value of y at t=0, lifetime in s
N, tau = 10000, 28
# Maximum time to consider (s)
tmax = 100
# A suitable grid of time points, and the exponential decay itself
t = np.linspace(0, tmax, 1000)
y = N * np.exp(-t/tau)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, y)
# The number of lifetimes that fall within the plotted time interval
ntau = tmax // tau + 1
# xticks at 0, tau, 2*tau, ..., ntau*tau; yticks at the corresponding y-values
xticks = [i*tau for i in range(ntau)]
yticks = [N * np.exp(-i) for i in range(ntau)]
ax.set_xticks(xticks)
ax.set_yticks(yticks)
# xtick labels: 0, tau, 2tau, ...
xtick_labels = [r'$0#39;, r'$\tau#39;] + [r'${}\tau#39;.format(k) for k in range(2,ntau)]
ax.set_xticklabels(xtick_labels)
# corresponding ytick labels: N, N/e, N/2e, ...
ytick_labels = [r'$N#39;,r'$N/e#39;] + [r'$N/{}e#39;.format(k) for k in range(2,ntau)]
ax.set_yticklabels(ytick_labels)
ax.set_xlabel(r'$t\;/\mathrm{s}#39;)
ax.set_ylabel(r'$y#39;)
ax.grid()
plt.show()
x轴标记被设置为1，τ，2τ，…y轴标记为N,N/e,N/2e，…
注意，标记序列的长度必须与所需的标记值列表的长度相对应。
注释图
另一个注释图的例子，这次是英国石油公司(BP plc, LSE: BP)的股价，其中添加了一些值得注意的事件。
数据来源 https://consent.yahoo.com/v2/collectConsent?sessionId=4_cc-session_adac4516-5607-4dd8-a394-f9ca47541e5b
import datetime
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.dates import strpdate2num
from datetime import datetime
def date_to_int(s):
epoch = datetime(year=1970, month=1, day=1)
date = datetime.strptime(s, '%Y-%m-%d')
return (date - epoch).days
def bindate_to_int(bs):
return date_to_int(bs.decode('ascii'))
dt = np.dtype([('daynum','i8'), ('close', 'f8')])
share_price = np.loadtxt('bp-share-prices.csv', skiprows=1, delimiter=',',
usecols=(0,4), converters={0: bindate_to_int},
dtype=dt)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(share_price['daynum'], share_price['close'], c='g')
ax.fill_between(share_price['daynum'], 0, share_price['close'], facecolor='g',
alpha=0.5)
daymin, daymax = share_price['daynum'].min(), share_price['daynum'].max()
ax.set_xlim(daymin, daymax)
price_max = share_price['close'].max()
def get_xy(date):
""" Return the (x,y) coordinates of the share price on a given date. """
x = date_to_int(date)
return share_price[np.where(share_price['daynum']==x)][0]
# A horizontal arrow and label
x,y = get_xy('1999-10-01')
ax.annotate('Share split', (x,y), xytext = (x+1000,y), va='center',
arrowprops=dict(facecolor='black', shrink=0.05))
# A vertical arrow and label
x,y = get_xy('2010-04-20')
ax.annotate('Deepwater Horizon\noil spill', (x,y), xytext = (x,price_max*0.9),
arrowprops=dict(facecolor='black', shrink=0.05), ha='center')
years = range(1989,2015,2)
ax.set_xticks([date_to_int('{:4d}-01-01'.format(year)) for year in years])
ax.set_xticklabels(years, rotation=90)
plt.show()
设置置信区域
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
FEMALE, MALE = 0, 1
dt = np.dtype([('mass', 'f8'), ('height', 'f8'), ('gender', 'i2')])
data = np.loadtxt('body.dat.txt', usecols=(22,23,24), dtype=dt)
fig, ax = plt.subplots()
def get_cov_ellipse(cov, centre, nstd, **kwargs):
"""
Return a matplotlib Ellipse patch representing the covariance matrix
cov centred at centre and scaled by the factor nstd.
"""
# Find and sort eigenvalues and eigenvectors into descending order
eigvals, eigvecs = np.linalg.eigh(cov)
order = eigvals.argsort()[::-1]
eigvals, eigvecs = eigvals[order], eigvecs[:, order]
# The anti-clockwise angle to rotate our ellipse by
vx, vy = eigvecs[:,0][0], eigvecs[:,0][1]
theta = np.arctan2(vy, vx)
# Width and height of ellipse to draw
width, height = 2 * nstd * np.sqrt(eigvals)
return Ellipse(xy=centre, width=width, height=height,
angle=np.degrees(theta), **kwargs)
labels, colours =['Female', 'Male'], ['magenta', 'blue']
for gender in (FEMALE, MALE):
sdata = data[data['gender']==gender]
height_mean = np.mean(sdata['height'])
mass_mean = np.mean(sdata['mass'])
cov = np.cov(sdata['mass'], sdata['height'])
ax.scatter(sdata['height'], sdata['mass'], color=colours[gender],
label=labels[gender], s=3)
e = get_cov_ellipse(cov, (height_mean, mass_mean), 3,
fc=colours[gender], alpha=0.4)
ax.add_artist(e)
ax.set_xlim(140, 210)
ax.set_ylim(30, 120)
ax.set_xlabel('Height /cm')
ax.set_ylabel('Mass /kg')
ax.legend(loc='upper left', scatterpoints=1)
plt.show()
该函数np.arctan2返回“双参数反正切”：np.arctan2(y, x)是正切之间的弧度角x-轴和点 ( x , y).
刻度线自定义
下面的程序创建一个带有主要和次要刻度线的绘图，自定义为比默认值更粗和更宽，主要刻度线指向绘图区域的内外。
import numpy as np
import matplotlib.pyplot as plt
# A selection of functions on rn abcissa points for 0 <= x < 1
rn = 100
rx = np.linspace(0, 1, rn, endpoint=False)
def tophat(rx):
""" Top hat function: y = 1 for x < 0.5, y=0 for x >= 0.5 """
ry = np.ones(rn)
ry[rx>=0.5]=0
return ry
# A dictionary of functions to choose from
ry = {'half-sawtooth': lambda rx: rx.copy(),
'top-hat': tophat,
'sawtooth': lambda rx: 2 * np.abs(rx-0.5)}
# Repeat the chosen function nrep times
nrep = 4
x = np.linspace(0, nrep, nrep*rn, endpoint=False)
y = np.tile(ry['top-hat'](rx), nrep)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(x,y, 'k', lw=2)
# Add a bit of padding around the plotted line to aid visualization
ax.set_ylim(-0.1,1.1)
ax.set_xlim(x[0]-0.5, x[-1]+0.5)
# Customize the tick marks and turn the grid on
ax.minorticks_on()
ax.tick_params(which='major', length=10, width=2, direction='inout')
ax.tick_params(which='minor', length=5, width=2, direction='in')
ax.grid(which='both')
plt.show()
该np.tile方法通过重复给定的数组nrep次数来构造一个数组。要绘制不同的周期函数，请在此处选择'half-sawtooth'或'sawtooth'。
饼图
温室气体排放
以下程序按“碳当量”的质量描述了温室气体的排放量。数据来自 2007 IPCC 报告。
import numpy as np
import matplotlib.pyplot as plt
# Annual greenhouse gas emissions, billion tons carbon equivalent (GtCe)
gas_emissions = np.array([(r'$\mathrm{CO_2}$-d', 2.2),
(r'$\mathrm{CO_2}$-f', 8.0),
('Nitrous\nOxide', 1.0),
('Methane', 2.3),
('Halocarbons', 0.1)],
dtype=[('source', 'U17'), ('emission', 'f4')])
# 5 colours beige
colours = ['#C7B299', '#A67C52', '#C69C6E', '#754C24', '#534741']
explode = [0, 0, 0.1, 0, 0]
fig, ax = plt.subplots()
ax.axis('equal') # So our pie looks round!
ax.pie(gas_emissions['emission'], colors=colours, shadow=True, startangle=90,
explode=explode, labels=gas_emissions['source'], autopct='%.1f%%',
pctdistance=1.15, labeldistance=1.3)
plt.show()
百分比值格式化为一位小数 ( autopct='%.1f%%')
分类柱状图
德国的可再生能源
该文件germany-energy-sources.txt包含 1990 年至 2018 年德国生产的可再生电力数据：
Renewable electricity generation in Germany in GWh (million kWh)
Year Hydro Wind Biomass Photovoltaics
2018 17974 109951 50851 45784
2017 20150 105693 50917 39401
2016 20546 79924 50928 38098
2015 18977 80624 50326 38726
2014 19587 58497 48287 36056
2013 22998 52737 45513 31010
2012 21755 51680 43203 26380
2011 17671 49857 36891 19599
2010 20953 38547 33924 11729
2009 19031 39420 30886 6583
2008 20443 41385 28014 4420
2007 21170 40507 24616 3075
2006 20031 31324 18934 2220
2005 19638 27774 14706 1282
2004 20745 26019 10636 557
2003 18322 19087 8948 313
2002 23124 16102 6048 162
2001 22733 10719 5214 76
2000 21732 9703 4731 60
1999 19647 5639 3585 30
1998 17216 4579 3256 35
1997 17357 3025 2273 18
1996 21957 2073 2098 12
1995 21780 1530 2010 7
1994 19930 927 1875 7
1993 17878 612 1635 3
1992 17397 281 1558 4
1991 14891 102 1471 1
1990 17426 72 1435 1Renewable electricity generation in Germany in GWh (million kWh)
Year Hydro Wind Biomass Photovoltaics
2013 21200 49800 47800 29300
2012 21793 50670 43350 26380
2011 17671 48883 37603 19559
...
下面的程序将这些数据绘制为堆积条形图，使用 Matplotlib 的填充图案来区分不同的来源。
import numpy as np
import matplotlib.pyplot as plt
data = np.loadtxt('germany-energy-sources.txt', skiprows=2, dtype='f8')
years = data[:,0]
n = len(years)
# GWh to TWh
data[:,1:] /= 1000
fig = plt.figure()
ax = fig.add_subplot(111)
sources = ('Hydroelectric', 'Wind', 'Biomass', 'Photovoltaics')
hatch = ['oo', '', 'xxxx', '//']
bottom = np.zeros(n)
bars = [None]*n
for i, source in enumerate(sources):
bars[i] = ax.bar(years, bottom=bottom, height=data[:,i+1], color='w',
hatch=hatch[i], align='center', edgecolor='k')
bottom += data[:,i+1]
ax.set_xticks(years)
plt.xticks(rotation=90)
ax.set_xlim(1989, 2019)
ax.set_ylabel('Renewable Electricity (TWh)')
ax.set_title('Renewable Electricity Generation in Germany, 1990-2018')
plt.legend(bars, sources, loc='best')
plt.show()
要包含图例，每个条形图对象必须存储在一个列表bars中，该列表通过ax.legend相应的标签序列传递给方法sources。
子图
该代码生成了10个子图，描述了n=0,1，⋯⋯9时sin(nπx)的图形。子地块间距的配置使它们在垂直方向上相互“碰撞”。
import numpy as np
import matplotlib.pyplot as plt
nrows = 10
fig, axes = plt.subplots(nrows,1)
# Zero vertical space between subplots
fig.subplots_adjust(hspace=0)
x = np.linspace(0,1,1000)
for i in range(nrows):
# n=nrows for the top subplot, n=0 for the bottom subplot
n = nrows - i
axes[i].plot(x, np.sin(n * np.pi * x), 'k', lw=2)
# We only want ticks on the bottom of each subplot
axes[i].xaxis.set_ticks_position('bottom')
if i < nrows-1:
# Set ticks at the nodes (zeros) of our sine functions
axes[i].set_xticks(np.arange(0, 1, 1/n))
# We only want labels on the bottom subplot xaxis
axes[i].set_xticklabels('')
axes[i].set_yticklabels('')
plt.show()
电磁频谱的描述
下面的程序使用text、axvline、axhline和axvspan注释一个简单的波图，以指示电磁频谱的不同区域。
import numpy as np
import matplotlib.pyplot as plt
# wavelength range, nm
lmin, lmax = 250, 1000
x = np.linspace(lmin, lmax, 1000)
# A wave with a smoothly increasing wavelength
wv = (np.sin(10 * np.pi * x / (lmax+lmin-x)))[::-1]
fig = plt.figure()
ax = fig.add_subplot(111, facecolor='k')
ax.plot(x, wv, c='w', lw=2)
ax.set_xlim(250,1000)
ax.set_ylim(-2,2)
# Label and delimit the different regions of the electromagnetic spectrum
ax.text(310, 1.5, 'UV', color='w', fontdict={'fontsize': 20})
ax.text(530, 1.5, 'Visible', color='k', fontdict={'fontsize': 20})
ax.annotate('', (400, 1.3), (750, 1.3), arrowprops={'arrowstyle': '<|-|>',
'color': 'w', 'lw': 2})
ax.text(860, 1.5, 'IR', color='w', fontdict={'fontsize': 20})
ax.axvline(400, -2, 2, c='w', ls='--')
ax.axvline(750, -2, 2, c='w', ls='--')
# Horizontal "axis" across the centre of the wave
ax.axhline(c='w')
# Ditch the y-axis ticks and labels; label the x-axis
ax.yaxis.set_visible(False)
ax.set_xlabel(r'$\lambda\;/\mathrm{nm}#39;)
# Finally, add some colourful rectangles representing a rainbow in the
# visible region of the spectrum.
# Dictionary mapping of wavelength regions (nm) to approximate RGB values
rainbow_rgb = { (400, 440): '#8b00ff', (440, 460): '#4b0082',
(460, 500): '#0000ff', (500, 570): '#00ff00',
(570, 590): '#ffff00', (590, 620): '#ff7f00',
(620, 750): '#ff0000'}
for wv_range, rgb in rainbow_rgb.items():
ax.axvspan(*wv_range, color=rgb, ec='none', alpha=1)
plt.show()
五颜六色的几何图形
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon, Circle, Rectangle
red, blue, yellow, green = '#ff0000', '#0000ff', '#ffff00', '#00ff00'
square = Rectangle((0.7, 0.1), 0.25, 0.25, facecolor=red)
circle = Circle((0.8, 0.8), 0.15, facecolor=blue)
triangle = Polygon(((0.05,0.1), (0.396,0.1), (0.223, 0.38)), fc=yellow)
rhombus = Polygon(((0.5,0.2), (0.7,0.525), (0.5,0.85), (0.3,0.525)), fc=green)
fig = plt.figure()
ax = fig.add_subplot(111, facecolor='k', aspect='equal')
for shape in (square, circle, triangle, rhombus):
ax.add_artist(shape)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
plt.show()
电偶极子的静电势
下面的代码给出了p=(qd,0,0)的电偶极子在(x,y)平面的电势等势面图，其中q=1.602×10−19C,d=1pm，使用点偶极近似。
import numpy as np
import matplotlib.pyplot as plt
# Dipole charge (C), Permittivity of free space (F.m-1)
q, eps0 = 1.602e-19, 8.854e-12
# Dipole +q, -q distance (m) and a convenient combination of parameters
d = 1.e-12
k = 1/4/np.pi/eps0 * q * d
# Cartesian axis system with origin at the dipole (m)
X = np.linspace(-5e-11, 5e-11, 1000)
Y = X.copy()
X, Y = np.meshgrid(X, Y)
# Dipole electrostatic potential (V), using point dipole approximation
Phi = k * X / np.hypot(X, Y)**3
fig = plt.figure()
ax = fig.add_subplot(111)
# Draw contours at values of Phi given by levels
levels = np.array([10**pw for pw in np.linspace(0,5,20)])
levels = sorted(list(-levels) + list(levels))
# Monochrome plot of potential
ax.contour(X, Y, Phi, levels=levels, colors='k', linewidths=2)
plt.show()
简单的等高线图
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
X = np.linspace(0,1,100)
Y = X.copy()
X, Y = np.meshgrid(X, Y)
alpha = np.radians(25)
cX, cY = 0.5, 0.5
sigX, sigY = 0.2, 0.3
rX = np.cos(alpha) * (X-cX) - np.sin(alpha) * (Y-cY) + cX
rY = np.sin(alpha) * (X-cX) + np.cos(alpha) * (Y-cY) + cY
Z = (rX-cX)*np.exp(-((rX-cX)/sigX)**2) * np.exp(- ((rY-cY)/sigY)**2)
fig = plt.figure()
ax = fig.add_subplot(111)
# Reversed Greys colourmap for filled contours
cpf = ax.contourf(X,Y,Z, 20, cmap=cm.Greys_r)
# Set the colours of the contours and labels so they're white where the
# contour fill is dark (Z < 0) and black where it's light (Z >= 0)
colours = ['w' if level<0 else 'k' for level in cpf.levels]
cp = ax.contour(X, Y, Z, 20, colors=colours)
ax.clabel(cp, fontsize=12, colors=colours)
plt.show()
分形图
巴恩斯利蕨类是一种分形，类似于蕨类中的黑色脾苔种。它是通过在(x,y)平面上绘制从(0,0)开始的一系列点来构建的，这些点由以下的仿射变换f1, f2, f3和f4生成，其中每个变换应用于前一个点，并随机选择概率p1=0.01, p2=0.85, p3=0.07和p4=0.07。.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
f1 = lambda x,y: (0., 0.16*y)
f2 = lambda x,y: (0.85*x + 0.04*y, -0.04*x + 0.85*y + 1.6)
f3 = lambda x,y: (0.2*x - 0.26*y, 0.23*x + 0.22*y + 1.6)
f4 = lambda x,y: (-0.15*x + 0.28*y, 0.26*x + 0.24*y + 0.44)
fs = [f1, f2, f3, f4]
npts = 50000
# Canvas size (pixels)
width, height = 300, 300
aimg = np.zeros((width, height))
x, y = 0, 0
for i in range(npts):
# Pick a random transformation and apply it
f = np.random.choice(fs, p=[0.01, 0.85, 0.07, 0.07])
x, y = f(x,y)
# Map (x,y) to pixel coordinates.
# NB we "know" that -2.2 < x < 2.7 and 0 <= y < 10
ix, iy = int(width / 2 + x * width / 10), int(y * height / 12)
# Set this point of the array to 1 to mark a point in the fern
aimg[iy, ix] = 1
plt.imshow(aimg[::-1,:], cmap=cm.Greens)
plt.show()
三维图——螺旋线
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
n = 1000
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# Plot a helix along the x-axis
theta_max = 8 * np.pi
theta = np.linspace(0, theta_max, n)
x = theta
z = np.sin(theta)
y = np.cos(theta)
ax.plot(x, y, z, 'b', lw=2)
# An line through the centre of the helix
ax.plot((-theta_max*0.2, theta_max * 1.2), (0,0), (0,0), color='k', lw=2)
# sin/cos components of the helix (e.g. electric and magnetic field
# components of a circularly-polarized electromagnetic wave
ax.plot(x, y, 0, color='r', lw=1, alpha=0.5)
ax.plot(x, [0]*n, z, color='m', lw=1, alpha=0.5)
# Remove axis planes, ticks and labels
ax.set_axis_off()
plt.show()
简单曲面图
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
L, n = 2, 400
x = np.linspace(-L, L, n)
y = x.copy()
X, Y = np.meshgrid(x, y)
Z = np.exp(-(X**2 + Y**2))
fig, ax = plt.subplots(nrows=2, ncols=2, subplot_kw={'projection': '3d'})
ax[0,0].plot_wireframe(X, Y, Z, rstride=40, cstride=40)
ax[0,1].plot_surface(X, Y, Z, rstride=40, cstride=40, color='m')
ax[1,0].plot_surface(X, Y, Z, rstride=12, cstride=12, color='m')
ax[1,1].plot_surface(X, Y, Z, rstride=20, cstride=20, cmap=cm.hot)
for axes in ax.flatten():
axes.set_xticks([-2, -1, 0, 1, 2])
axes.set_yticks([-2, -1, 0, 1, 2])
axes.set_zticks([0, 0.5, 1])
fig.tight_layout()
plt.show()
动画
弹跳球动画
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Acceleration due to gravity, m.s-2.
g = 9.81
# The maximum x-range of ball's trajectory to plot.
XMAX = 5
# The coefficient of restitution for bounces (-v_up/v_down).
cor = 0.65
# The time step for the animation.
dt = 0.005
# Initial position and velocity vectors.
x0, y0 = 0, 4
vx0, vy0 = 1, 0
def get_pos(t=0):
"""A generator yielding the ball's position at time t."""
x, y, vx, vy = x0, y0, vx0, vy0
while x < XMAX:
t += dt
x += vx0 * dt
y += vy * dt
vy -= g * dt
if y < 0:
# bounce!
y = 0
vy = -vy * cor
yield x, y
def init():
"""Initialize the animation figure."""
ax.set_xlim(0, XMAX)
ax.set_ylim(0, y0)
ax.set_xlabel('$x$ /m')
ax.set_ylabel('$y$ /m')
line.set_data(xdata, ydata)
ball.set_center((x0, y0))
height_text.set_text(f'Height: {y0:.1f} m')
return line, ball, height_text
def animate(pos):
"""For each frame, advance the animation to the new position, pos."""
x, y = pos
xdata.append(x)
ydata.append(y)
line.set_data(xdata, ydata)
ball.set_center((x, y))
height_text.set_text(f'Height: {y:.1f} m')
return line, ball, height_text
# Set up a new Figure, with equal aspect ratio so the ball appears round.
fig, ax = plt.subplots()
ax.set_aspect('equal')
# These are the objects we need to keep track of.
line, = ax.plot([], [], lw=2)
ball = plt.Circle((x0, y0), 0.08)
height_text = ax.text(XMAX*0.5, y0*0.8, f'Height: {y0:.1f} m')
ax.add_patch(ball)
xdata, ydata = [], []
interval = 1000*dt
ani = animation.FuncAnimation(fig, animate, get_pos, blit=True,
interval=interval, repeat=False, init_func=init)
plt.show()衰减正弦曲线
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
# Time step for the animation (s), max time to animate for (s).
dt, tmax = 0.01, 5
# Signal frequency (s-1), decay constant (s-1).
f, alpha = 2.5, 1
# These lists will hold the data to plot.
t, M = [], []
# Draw an empty plot, but preset the plot x- and y-limits.
fig, ax = plt.subplots()
line, = ax.plot([], [])
ax.set_xlim(0, tmax)
ax.set_ylim(-1, 1)
ax.set_xlabel('t /s')
ax.set_ylabel('M (arb. units)')
def init():
return line,
def animate(i, t, M):
"""Draw the frame i of the animation."""
# Append this time point and its data and set the plotted line data.
_t = i*dt
t.append(_t)
M.append(np.sin(2*np.pi*f*_t) * np.exp(-alpha*_t))
line.set_data(t, M)
return line,
# Interval between frames in ms, total number of frames to use.
interval, nframes = 1000 * dt, int(tmax / dt)
# Animate once (set repeat=False so the animation doesn't loop).
ani = animation.FuncAnimation(fig, animate, frames=nframes, init_func=init,
fargs=(t, M), repeat=False, interval=interval, blit=True)
plt.show()
注意：分配给fargs参数的任何对象FuncAnimation都将传递给动画函数。
热力图
import numpy as np
import matplotlib.pyplot as plt
# Read in the relevant data from our input file.
dt = np.dtype([('month', np.int), ('day', np.int), ('T', np.float)])
data = np.genfromtxt('boston2019.dat', dtype=dt, usecols=(1, 2, 3),
delimiter=(4, 2, 2, 6))
# In our heatmap, nan will mean "no such date", e.g. 31 June.
heatmap = np.empty((12, 31))
heatmap[:] = np.nan
for month, day, T in data:
# NumPy arrays are zero-indexed; days and months are not!
heatmap[month-1, day-1] = T
# Plot the heatmap, customize and label the ticks.
fig = plt.figure()
ax = fig.add_subplot()
im = ax.imshow(heatmap, interpolation='nearest')
ax.set_yticks(range(12))
ax.set_yticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
days = np.array(range(0, 31, 2))
ax.set_xticks(days)
ax.set_xticklabels(['{:d}'.format(day+1) for day in days])
ax.set_xlabel('Day of month')
ax.set_title('Maximum daily temperatures in Boston, 2019')
# Add a colorbar along the bottom and label it.
cbar = fig.colorbar(ax=ax, mappable=im, orientation='horizontal')
cbar.set_label('Temperature, $^\circ\mathrm{C}#39;)
plt.show()
传递给的“可映射”对象fig.colorbar是由AxesImage返回的对象ax.imshow。
声明:本文由"麦兜"发布，不代表"破站"立场，转载联系作者并注明出处：https://www.848586.cn/137826.html
