数学曲线

一、蝴蝶曲线

$\begin{aligned} x=\sin(t)(e^{\cos(t)}-2\cos(4t)-\sin^5(t/12)) \\ y=\cos(t)(e^{\cos(t)}-2\cos(4t)-\sin^5(t/12)) \end{aligned},0\le t\le 12\pi$

极坐标:

$r=e^{\sin(\theta)}-2\cos(4\theta)+\sin(\frac{2\theta-\pi}{24})$

python代码:

import numpy as np
import matplotlib.pyplot as plt

t = np.arange(0.0, 12*np.pi, 0.01)
x = np.sin(t)*(np.e**np.cos(t) - 2*np.cos(4*t)-np.sin(t/12)**5)
y = np.cos(t)*(np.e**np.cos(t) - 2*np.cos(4*t)-np.sin(t/12)**5)

plt.figure(figsize=(8,6))
plt.axis('off')
plt.plot(x,y,color='blue',linewidth = '2')
plt.savefig("butter.jpg",dpi=400)

二、心形曲线

$(x^2+y^2-1)^3-x^2y^3=0$

python程序:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.font_manager import FontProperties
font=FontProperties(fname='/Library/Fonts/Songti.ttc',size=10)

x_coords = np.linspace(-100, 100, 500)
y_coords = np.linspace(-100, 100, 500)
points = [ ]
for y in y_coords:
    for x in x_coords:
        if ((x*0.03)**2+(y*0.03)**2-1)**3-(x*0.03)**2*(y*0.03)**3 <= 0:
            points.append({"x": x, "y": y})
heart_x = list(map(lambda point: point["x"], points))
heart_y = list(map(lambda point: point["y"], points))


plt.figure("heart", figsize=(6,15))
plt.subplot(521)
plt.scatter(heart_x, heart_y, s=10, alpha=0.5)
plt.title("普通",fontproperties=font)
plt.axis('off')

cmap =[("autumn","橙色-热情洋溢"),("cool","紫色-优雅宁静"),("magma","晚霞-醇厚脱俗"),("rainbow","彩虹-绚丽幻想"),("Reds","炽热-热烈奔放"),("spring","青春-充满朝气"),("viridis","翡翠-平静柔和"),("gist_rainbow","五彩缤纷-多姿多彩")]

for i in range(len(cmap)):
    plt.subplot(522+i)
    plt.scatter(heart_x, heart_y, s=10, alpha=0.5, c=range(len(heart_x)), cmap=cmap[i][0])
    plt.title(cmap[i][1],fontproperties=font)
    plt.axis('off')
# plt.show()
plt.savefig("heart.jpg",dpi=400)

三、Mandelbrot集合

曼德博集合可以用复二次多项式来定义:

$f_{c}(z)=z^{2}+c$

其中 $c$ 是一个复数参数。

从 $z=0$ 开始对 $f_{c}(z)$ 进行迭代:

$z_{n+1}=z_{n}^{2}+c,n=0,1,2,...$

$z_{0}=0$

$z_{1}=z_{0}^{2}+c=c$

$z_{2}=z_{1}^{2}+c=c^{2}+c$

每次迭代的值依序如以下序列所示: $(0,f_{c}(0),f_{c}(f_{c}(0)),f_{c}(f_{c}(f_{c}(0))),\ldots )$

不同的参数 $c$ 可能使序列的绝对值逐渐发散到无限大,也可能收敛在有限的区域内。曼德博集合 $M$ 就是使序列不延伸至无限大的所有复数 $c$ 的集合。

曼德博集合有一些特性和定理,可参考:曼德博集合

python代码

import numpy as np
import matplotlib.pylab as plt

MAXITERS,RADIUS = 200,100

def color(z, i):
    v = np.log2(i + 1 - np.log2(np.log2(abs(z)))) / 5
    if v < 1.0: return v**4, v**2.5, v
    else:
        v = max(0, 2-v)
        return v, v**1.5, v**3

def iterate(c):
    z = 0j
    for i in range(MAXITERS):
        if z.real*z.real + z.imag*z.imag > RADIUS: return color(z, i)
        z = z*z + c
    return 0, 0, 0

xmin,xmax, ymin,ymax,width,height=-2.1, 0.8, -1.16, 1.16, 800, 640
y, x = np.ogrid[ymin: ymax: height*1j, xmin: xmax: width*1j]
z = x + y*1j
R, G, B = np.asarray(np.frompyfunc(iterate, 1, 3)(z)).astype(np.float)
img = np.uint8(np.stack((R, G, B), axis=2) * 255)
plt.figure("Mandelbrot")
plt.imshow(img)
plt.axis('off')
plt.title('Mandelbrot')
# plt.show()
plt.savefig("Mandelbrot.jpg",dpi=400)

四、Newton分形

牛顿分形就是分形的一种,它与解方程的牛顿法(跟优化中的牛顿法是不同的方法)紧密相关。假如需要求解方程 $f(x) = 0$ 。

其中,x的定义域是整个复平面。如何求解这个方程的解呢？我们可以用牛顿法。牛顿法是一种数值解法,我们首先会估算一个“比较好”的初始值 $x_0$ ,然后使用迭代公式:

$x_{n+1} = x_n - \frac{f(x_n)}{f^{'}(x_n)}$

牛顿法可以确保,如果初始猜测值在根附近,那么迭代必然收敛。而且牛顿法是个二阶方法,收敛速度相当的快。下图是迭代一步的示意图

在 $x_1$ 处沿着方向 $\frac{f(x_1)}{f'(x_1)}$ 下降,与 $x$ 轴的交点即为 $x_2$ ,循环往复就能得到方程的根。

学过中学数学的我们都知道, n次方程在复数域上有n个根,那么用牛顿法收敛的根就可能有n个目标。牛顿法收敛到哪个根取决于迭代的起始值。根据最后的收敛结果,我们把所有收敛到同一个根的起始点画上同一种颜色,最终就形成了牛顿分形图。

python代码

import numpy as np
import matplotlib.pyplot as plt
from numba import jit

# define functions manually, do not use numpy's poly1d funciton!
@jit('complex64(complex64)', nopython=True)
def f(z):
    # z*z*z is faster than z**3
    # return z*z*z - 1
    return z**5 + 0.25*z*z + 1.17

@jit('complex64(complex64)', nopython=True)
def df(z):
    # return 3*z*z
    return 5*z**4 + 0.5*z

@jit('float64(complex64)', nopython=True)
def iterate(z):
    num = 0
    while abs(f(z)) > 1e-4:
        w = z - f(z)/df(z)
        num += np.exp(-1/abs(w-z))
        z = w
    return num

imgsize = 600
y, x = np.ogrid[1: -1: imgsize*2j, -1: 1: imgsize*2j]
z = x + y * 1j
img = np.frompyfunc(iterate, 1, 1)(z).astype(np.float)
fig = plt.figure(figsize=(imgsize/100.0, imgsize/100.0), dpi=100)
ax = fig.add_axes([0, 0, 1, 1], aspect=1)
ax.axis('off')
ax.imshow(img, cmap='hot')
fig.savefig('newton.jpg')

五、洛伦茨吸引子

含时间参数的形式:

$\begin{cases} & \frac{dx(t)}{dt}=\sigma(y(t)-x(t)) \\ & \frac{dy(t)}{dt}=\rho x(t)-y(t)-x(t)z(t) \\ & \frac{dz(t)}{dt}=x(t)y(t)-\beta z(t) \end{cases}$

$\sigma$ 称为普兰特尔数, $\rho$ 称为瑞利数。所有的 $\sigma,\rho,\beta>0$ ,但通常 $\sigma=10,\beta=8/3,\rho$ 不定。若 $\rho<1$ ,则吸引子为原点,没有任何其他稳定点。 $1≤\rho<13.927$ 时,螺线轨迹接近两点(这相当于存在阻尼振子),两点的位置由下列式子决定: $x=\pm \sqrt{b(\rho-1)},y=\pm \sqrt{b(\rho-1)}、z=\rho-1$ 。

python代码:

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
from scipy.integrate import odeint
from mpl_toolkits.mplot3d import Axes3D

# number of particles.
num_particles = 20

# constants for Lorenz system.
alpha, beta,gamma = 10.0, 8/3.0, 28.0

def derivative(point, t):
    # return the tangent direction at (x,y,z)
    x, y, z = point
    return [alpha * (y - x),x * (gamma - z) - y,x * y - beta * z]

fig = plt.figure(figsize=(6.4, 4.8), dpi=100)
ax = fig.add_axes([0, 0, 1, 1], projection='3d', xlim=(-25, 25), ylim=(-35, 35), zlim=(5, 55))
ax.view_init(30, 0)
ax.axis('off')

lines, points = [], []
colors = plt.cm.gist_ncar(np.linspace(0, 1, num_particles))

for c in colors:
    lines.extend(ax.plot([], [], '-', c=c))
    points.extend(ax.plot([], [], 'o', c=c))

x0 = -15 + 30 * np.random.random((num_particles, 3))
t = np.linspace(0, 40, 10001)
x_t = np.array([odeint(derivative, point, t) for point in x0])

def init():
    for line, point in zip(lines, points):
        line.set_data([], [])
        line.set_3d_properties([])

        point.set_data([], [])
        point.set_3d_properties([])
    return lines + points

def animate(i):
    i = 2*i % x_t.shape[1]  # accelarate the animation.

    for line, point, x_j in zip(lines, points, x_t):
        x, y, z = x_j[:i].T
        line.set_data(x, y)
        line.set_3d_properties(z)

        # note that plot() receives a list parameter so we have
        # to write x[-1:] instead of x[-1]!
        point.set_data(x[-1:], y[-1:])
        point.set_3d_properties(z[-1:])

    ax.view_init(30, 0.3*i)
    fig.canvas.draw()
    return lines + points


anim = FuncAnimation(fig, animate, init_func=init, interval=20, frames=2000, blit=True)
anim.save('lorenz.mp4', writer='ffmpeg', fps=30, dpi=200,bitrate=1000, codec='libx264', extra_args=['-crf', '10'])

六、羊角螺线

羊角螺线(clothoid),又称欧拉螺线(Euler spiral),是形式为

$x=C(t)$

$y=S(t)$

的曲线,其中 $C(t)$ 、 $S(t)$ 为Fresnel积分:

$S(x)=\int _{0}^{x}\sin(t^{2})\,dt=\sum _{{n=0}}^{{\infty }}(-1)^{n}{\frac {x^{{4n+3}}}{(4n+3)(2n+1)!}},$

$C(x)=\int _{0}^{x}\cos(t^{2})\,dt=\sum _{{n=0}}^{{\infty }}(-1)^{n}{\frac {x^{{4n+1}}}{(4n+1)(2n)!}}.$

上面参数方程的参数 $t$ ,也是螺线于该点的曲率: $\kappa (t)=2t$ 。两个螺线的中心位于 $\frac {\sqrt {\pi }}{2}$ , $\frac {\sqrt {\pi }}{2}$

由于此螺线的曲率与长度成正比,故常用于公路工程或铁路工程,以缓和直路线与圆曲路线之间的曲率变化(向心力变化)。在光学上,近场衍射(Fresnel衍射)中会应用Fresnel积分。

python代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
import math

t = np.arange(-10,10,0.01)
s = map(lambda i:integrate.quad(lambda x:math.sin(pow(x,2)),0,i),t)
c = map(lambda i:integrate.quad(lambda x:math.cos(pow(x,2)),0,i),t)

s = [i[0] for i in list(s)]
c = [i[0] for i in list(c)]

plt.figure(figsize=(4,3))
plt.axis('off')
plt.plot(c,s, color='red',linewidth = '1')
plt.savefig("butter.jpg",bbox_inches='tight', transparent=True,dpi=400)