Python学数学（一）

用几何学变换形状

translate(50,80)
# 会将整个网格向右移动50像素，向下移动80像素

#将图形绕着原点进行旋转
rotate()    #参数为弧度
rotate(radians(20)) #将角度转换为弧度

画一个圆圈图案

def setup():
    size(600,600)
    
def draw():
    translate(width/2,height/2)
    fill(0)
    for i in range(12):
        ellipse(200,0,50,50)
        rotate(radians(360/12))

画一圈方形

def setup():
    size(600,600)
    
def draw():
    translate(width/2,height/2)
    fill(0)
    for i in range(12):
        rect(200,0,50,50)
        rotate(radians(360/12))

使用pushMatrix()和popMatrix()来保存当前的位置

​
t = 0
​
def setup():
    size(600,600)
    
def draw():
    global t
    
    background(255)
    translate(width/2,height/2)
    rotate(radians(t))
    
    rectMode(CENTER)    #设置绘制方形的中心点为中心  CORNER 右上角  CORNERS 左下角 RADIUS 左上角
    
    for i in range(12):
        pushMatrix()    #开始旋转正方形之前保存原来的中心位置
        translate(200,0)
        rotate(radians(t))
        rect(0,0,50,50)
        popMatrix()     #恢复之前的位置
        rotate(radians(360/12))
    t += 0.1

用三角学制造振荡

绘制sin函数的图像

xmin = -10
xmax = 10
​
ymin = -10
ymax = 10
​
rangex = xmax - xmin
rangey = ymax - ymin
​
def setup():
    global xscl,yscl
    size(600,600)
    
    xscl = width / rangex
    yscl = -height / rangey
    
def draw():
    global xscl,yscl
    background(255)
    translate(width/2,height/2)
    grid(xscl,yscl)
    
    stroke(255,0,0)
    x = xmin
    while x <= xmax:
        fill(0)
        line(x * xscl,f(x)*yscl,(x+0.1)*xscl,f(x+0.1)*yscl)
        x += 0.1
    
​
def grid(xscl,yscl):
    strokeWeight(1)
    stroke(0,255,255)
    
    for i in range(xmin,xmax + 1):
        line(i,ymin,i,ymax)
        line(i * xscl,ymin * yscl,i * xscl,ymax * yscl)
        line(xmin * xscl,i * yscl,xmax * xscl,i * yscl)
        
    stroke(0)
    line(0,ymax * yscl ,0,ymin * yscl)
    line(xmin * xscl,0,xmax * xscl,0)
    
def f(x):
    return sin(x)

绘制多边形

def setup():
    size(600,600)
​
def draw():
    beginShape()    #定义形状的开始
    vertex(100,100) #利用vertex函数，提供各个点的坐标
    vertex(100,200)
    vertex(200,200)
    vertex(200,100)
    vertex(150,50)
    endShape(CLOSE) #定义形状的结束

绘制任意多边形

def setup():
    size(600,600)
    
    
def draw():
    polygon(100,5)
​
​
def polygon(r,number):
    translate(width/2,height/2)
    beginShape()
    
    for i in range(number):
        vertex(r*cos(radians(360/number*i)),r * sin(radians(360/number * i)))   #利用极坐标的方式绘制图形
    endShape(CLOSE)

绘制正弦波

r1 = 50
r2 = 10
t = 0
circleList = []

def setup():
    size(600,600)
    
def draw():
    global t,circleList
    background(200)
    
    translate(width/4,height/2)	#绘制大圆
    noFill()
    stroke(0)
    ellipse(0,0,2 * r1,2 * r1)
    
    x = r1 * cos(radians(t))	#计算小圆的坐标
    y = r1 * sin(radians(t))
    
    fill(255,0,0)				#绘制小圆
    ellipse(x,y,2 * r2 ,2 * r2)
    
    stroke(0,255,0)				#绘制绿色小圆与大圆上小圆之间的绿色线
    line(x,y,150,y)
    fill(0,255,0)
    ellipse(150,y,5,5)
    
    circleList = [y] + circleList[:249]	#记录每次绿色小圆走过的y坐标
    
    for i in range(len(circleList)):	#遍历小圆的路径，绘制一系列绿色圆
        ellipse(150 + i,circleList[i],5,5)
    
    t += 3

利用函数来获取列表的索引和值

绘制万花尺程序

不断改变小圆的中心位置，让小圆贴着大圆的边缘进行旋转

#定义需要的参数
r1 = 300
r2 = 175
r3 = 5.0

x1 = 0
y1 = 0
t = 0  #时间变量

points = []#用于存放空列表

def setup():
    size(600,600)
    
def draw():
    global t
    translate(width/2,height/2)
    background(255)
    
    noFill()
    stroke(0)
    
    #利用极坐标计算小圆的圆心位置
    x2 = (r1 - r2) * cos(radians(t))
    y2 = (r1 - r2) * sin(radians(t))
    
    ellipse(x1,y1,r1*2,r1*2)  #绘制大圆
    ellipse(x2,y2,r2 *2,r2 *2) #绘制小圆
    
    t += 0.8
    

在小圆中心在添加一个小圆，用于绘制路径

#定义需要的参数
r1 = 300
r2 = 175
r3 = 5.0

x1 = 0
y1 = 0
t = 0  #时间变量

points = []#用于存放空列表

def setup():
    size(600,600)
    
def draw():
    global t
    translate(width/2,height/2)
    background(255)
    
    noFill()
    stroke(0)
    
    #利用极坐标计算小圆的圆心位置
    x2 = (r1 - r2) * cos(radians(t))
    y2 = (r1 - r2) * sin(radians(t))
    
    #计算画笔的位置
    x3 = x2 + (r2 - r3) * cos(radians(t))
    y3 = y2 + (r2 - r3) * sin(radians(t))
    
    ellipse(x1,y1,r1 * 2,r1 * 2)  #绘制大圆
    ellipse(x2,y2,r2 * 2,r2 * 2) #绘制小圆
    ellipse(x3,y3,r3 * 2,r3 * 2) #绘制画笔
    
    t += 0.8
    

记录点的位置，绘制直线

#定义需要的参数
r1 = 300
r2 = 105
r3 = 5.0
prop = 0.8

x1 = 0
y1 = 0
t = 0  #时间变量

points = []#用于存放空列表

def setup():
    size(600,600)
    
def draw():
    global t
    global points
    translate(width/2,height/2)
    background(255)
    
    noFill()
    stroke(0)
    
    #利用极坐标计算小圆的圆心位置
    x2 = (r1 - r2) * cos(radians(t))
    y2 = (r1 - r2) * sin(radians(t))
    
    #计算画笔的位置
    # x3 = x2 - prop * (r2 - r3) * cos(radians(-((r1 - r2) /r2) * t * 2))
    # y3 = y2 - prop * (r2 - r3) * sin(radians(-((r1 - r2) /r2) * t * 2))
    x3 =  prop * (r2 - r3) * cos(radians(t * 5)) + x2 
    y3 =  prop * (r2 - r3) * sin(radians(t * 5)) + y2
    
    points = [[x3,y3]] + points[:2000]
    

    ellipse(x1,y1,r1 * 2,r1 * 2)  #绘制大圆
    ellipse(x2,y2,r2 * 2,r2 * 2) #绘制小圆
    
    fill(255,0,0)
    ellipse(x3,y3,r3 * 2,r3 * 2) #绘制画笔
    
    #进行遍历，绘制直线
    for i,p in enumerate(points):
        if i < len(points) - 1:
            stroke(255,0,0)
            line(p[0],p[1],points[i + 1][0],points[i + 1][1])
    
    t += 0.8
    

复数

复数相加

def cAdd(a,b):
    return [a[0] + b[0],a[1] + b[1]]

print(cAdd([1,2],[3,4]))

复数相乘

(a + bi)(c + di) = ac - bd + (ad +bc)i

def cMult(u,v):
    return [u[0] * v[0]-u[1]*v[1],u[1]*v[0]+u[0]*v[1]]

复数的大小

复数的大小或者说是绝对值，是其在复数坐标平面上到原点的距离

import math as m

def magnitude(z):
    return m.sqrt(z[0]**2+z[1]**2)

print(magnitude([2,1]))

import math as m

def magnitude(z):
    return m.sqrt(z[0]**2+z[1]**2)

def cAdd(a,b):
    return [a[0] + b[0],a[1] + b[1]]

def cMult(u,v):
    return [u[0] * v[0]-u[1]*v[1],u[1]*v[0]+u[0]*v[1]]


def mandelrbot(z,num):
    count = 0

    z1 = z

    while count <= num:
        #如果超过了临界值，则返回迭代次数
        if magnitude(z1) > 2.0:
            return count

        z1 = cAdd(cMult(z1,z1),z)
        print(magnitude(z1))
        count += 1

    #如果一直没有超过临界值，则返回参数num
    return num

mandelrbot([0.25,0.75],5)

芒德布罗集

要绘制芒德布罗集，需要将这个画面的像素左边，转换成复数z，然后重复地将其平方加上原来的z。如果最后发现最后计算出来的z的值是一个发散的数据，则将这个点涂成白色，如果是一个收敛的值，则将它的这个点涂成黑色

import math as m
​
xmin = -2
xmax = 2
​
ymin = -2
ymax = 2
​
rangex = xmax - xmin
rangey = ymax - ymin
​
xscl = 0
yscl = 0
​
def magnitude(z):
    return m.sqrt(z[0]**2+z[1]**2)
​
def cAdd(a,b):
    return [a[0] + b[0],a[1] + b[1]]
​
def cMult(u,v):
    return [u[0] * v[0]-u[1]*v[1],u[1]*v[0]+u[0]*v[1]]
​
​
def mandelrbot(z,num):
    count = 0
​
    z1 = z
​
    while count <= num:
        #如果超过了临界值，则返回迭代次数
        if magnitude(z1) > 2.0:
            return count
        #求出新的复数：新复数的平方加上原复数
        z1 = cAdd(cMult(z1,z1),z)
        count += 1
    return num
​
​
def setup():
    global xscl,yscl
    size(600,600)
    
    #设置颜色模式为HSB
    colorMode(HSB)
    noStroke()
    xscl = float(rangex)/width
    yscl = float(rangey)/height
​
def draw():
    global points,xscl,yscl
    background(255)
    
    #遍历所有像素
    for x in range(width):
        for y in range(height):
            z = [(xmin + x * xscl),(ymin + y * yscl)]
            
            #看看这个像素对应的复数，是否发散
            col = mandelrbot(z,100)
            
            #如果不发散，就涂成黑色
            if col == 100:
                fill(0)
            
            #如果发散，就涂成白色
            else:
                fill(3 * col,255,255)
            rect(x,y,1,1)
    
​
​

茹利亚集

茹利亚集的构建方法：

求出复数的平方后不是加上像素对应的原复数，而是加上一个复数常数C，对于所有像素，C的值都是相同。通过使用不同的C值，我们就可以创建出许多茹利亚集

import math as m
​
xmin = -2
xmax = 2
​
ymin = -2
ymax = 2
​
rangex = xmax - xmin
rangey = ymax - ymin
​
xscl = 0
yscl = 0
​
def magnitude(z):
    return m.sqrt(z[0]**2+z[1]**2)
​
def cAdd(a,b):
    return [a[0] + b[0],a[1] + b[1]]
​
def cMult(u,v):
    return [u[0] * v[0]-u[1]*v[1],u[1]*v[0]+u[0]*v[1]]
​
​
def mandelrbot(z,c,num):
    count = 0
​
    z1 = z
​
    while count <= num:
        #如果超过了临界值，则返回迭代次数
        if magnitude(z1) > 2.0:
            return count
​
        z1 = cAdd(cMult(z1,z1),c)
        count += 1
    return num
​
​
def setup():
    global xscl,yscl
    size(600,600)
    
    #设置颜色模式为HSB
    colorMode(HSB)
    noStroke()
    xscl = float(rangex)/width
    yscl = float(rangey)/height
​
def draw():
    global points,xscl,yscl
    background(255)
    
    c = [0.285,0.01]
    
    #遍历所有像素
    for x in range(width):
        for y in range(height):
            z = [(xmin + x * xscl),(ymin + y * yscl)]
            
            #看看这个像素对应的复数，是否发散
            col = mandelrbot(z,c,100)
            
            #如果不发散，就涂成黑色
            if col == 100:
                fill(0)
            
            #如果发散，就涂成白色
            else:
                fill(3 * col,255,255)
            rect(x,y,1,1)
    
​
​

将矩阵用于计算机图形和方程组

矩阵相加

相同维度的矩阵可以进行相加运算，编写一个程序，来对相同维度的矩阵进行加法运算

def addMatrices(a,b):
    c = [[a[0][0] + b[0][0],a[0][1] + b[0][1]],
         [a[1][0] + b[1][0],a[1][1] + b[1][1]]]
    return c
    
A = [[2,3],[5,-8]]
B = [[1,-4],[8,-6]]
     
C = addMatrices(A,B)
print(C)

矩阵相乘

要求：第一个矩阵的行要等于第二个矩阵的列数

#a为取行，b为取列
def multmatrix(a,b):
    m = len(a)
    n = len(b[0])
    
    newmatrix = []
    
    for i in range(m):
        row = []
        for j in range(n):
            sum1= 0
            for k in range(len(b)):
                sum1 += a[i][k] * b[k][j]
            row.append(sum1)
        newmatrix.append(row)
    return newmatrix
​
B = [[1,2,-3,-1]]
A = [[4,-1],[-2,3],[6,-3],[1,0]]
​
print(multmatrix(B,A))

矩阵A*B与B*A是两个不同的结果，矩阵不支持交换律

变换矩阵

​
xmin = -10
xmax = 10
​
​
ymin = -10
ymax = 10
​
rangex = xmax - xmin
rangey = ymax - ymin
​
#一个字母F矩阵
fmatrix = [[0,0],[1,0],[1,2],[2,2],[2,3],[1,3],[1,4],[3,4],[3,5],[0,5],[0,0]]
​
def setup():
    global xscl, yscl
    size(600,600)
    xscl= width/rangex
    yscl= -height/rangey
    noFill()
​
def draw():
    global xscl, yscl,fmatrix
    background(255)
    translate(width/2,height/2)
    grid(xscl, yscl)
    
    strokeWeight(2)
    stroke(0)
    graphPoints(fmatrix)
    
    
#用线段连接列表中的相邻的点
#利用vertex函数，来连接每一个点
def graphPoints(matrix):
    beginShape()
    for pt in matrix:
        vertex(pt[0] * xscl,pt[1] * yscl)
    endShape()
​
def grid(xscl,yscl):
    strokeWeight(1)
    stroke(0,255,255)
    for i in range(xmin,xmax+1):
        line(i*xscl,ymin*yscl,i*xscl,ymax*yscl)
    for i in range(ymin,ymax+1):
        line(xmin*xscl,i*yscl,xmax*xscl,i*yscl)
    stroke(0)
    line(0,ymin*yscl,0,ymax*yscl)
    line(xmin*xscl,0,xmax*xscl,0)

绘制结果如下：

下面利用矩阵乘法，来对矩阵进行变换，我们需要一个变换矩阵，如下：$$
R_0 = \begin{bmatrix} cos\theta & -sin\theta \\ sin\theta & cos\theta \end{bmatrix}
$$

\theta为要旋转的角度，如果要逆时针旋转90\degree，由于cos(90)=0且sin(90)=1，那么获得的转换矩阵如下：$$
R = \begin{bmatrix} 0&-1\\1&0 \end{bmatrix}
$$

相乘之后，绘制出新的结果：

​
xmin = -10
xmax = 10
​
​
ymin = -10
ymax = 10
​
rangex = xmax - xmin
rangey = ymax - ymin
​
fmatrix = [[0,0],[1,0],[1,2],[2,2],[2,3],[1,3],[1,4],[3,4],[3,5],[0,5],[0,0]]
transformation_matrix = [[0,-1],[1,0]]
​
def setup():
    global xscl, yscl
    size(600,600)
    xscl= width/rangex
    yscl= -height/rangey
    noFill()
​
def draw():
    global xscl, yscl,fmatrix
    background(255)
    translate(width/2,height/2)
    grid(xscl, yscl)
    
    #绘制原来的矩阵
    strokeWeight(2)
    stroke(0)
    graphPoints(fmatrix)
    
    #绘制出新的变换后的矩阵
    strokeWeight(2)
    stroke(255,0,0)
    graphPoints(multmatrix(fmatrix,transformation_matrix))
    
    
#用线段连接列表中的相邻的点
def graphPoints(matrix):
    beginShape()
    for pt in matrix:
        vertex(pt[0] * xscl,pt[1] * yscl)
    endShape()
​
def grid(xscl,yscl):
    strokeWeight(1)
    stroke(0,255,255)
    for i in range(xmin,xmax+1):
        line(i*xscl,ymin*yscl,i*xscl,ymax*yscl)
    for i in range(ymin,ymax+1):
        line(xmin*xscl,i*yscl,xmax*xscl,i*yscl)
    stroke(0)
    line(0,ymin*yscl,0,ymax*yscl)
    line(xmin*xscl,0,xmax*xscl,0)
    
#a为取行，b为取列
def multmatrix(a,b):
    m = len(a)
    n = len(b[0])
    
    newmatrix = []
    
    for i in range(m):
        row = []
        for j in range(n):
            sum1= 0
            for k in range(len(b)):
                sum1 += a[i][k] * b[k][j]
            row.append(sum1)
        newmatrix.append(row)
    return newmatrix

但是发现它是一个顺时针旋转的过程，公认的做法是变换矩阵在前，而被转换的矩阵在后面$$
v^{‘} = R_\theta v_\theta
$$

因此，需要将F矩阵进行转置运算

#转置函数
def transpose(a):
    newmatrix = []
    for i in range(len(a[0])):
        row = []
        for j in range(len(a)):
            row.append(a[j][i])
        newmatrix.append(row)
    return newmatrix
​

转置后进行相乘，将结果在进行一次转置，得到正确结果

​
xmin = -10
xmax = 10
​
​
ymin = -10
ymax = 10
​
rangex = xmax - xmin
rangey = ymax - ymin
​
fmatrix = [[0,0],[1,0],[1,2],[2,2],[2,3],[1,3],[1,4],[3,4],[3,5],[0,5],[0,0]]
transformation_matrix = [[0,-1],[1,0]]
​
def setup():
    global xscl, yscl
    size(600,600)
    xscl= width/rangex
    yscl= -height/rangey
    noFill()
​
def draw():
    global xscl, yscl,fmatrix
    background(255)
    translate(width/2,height/2)
    grid(xscl, yscl)
    
    strokeWeight(2)
    stroke(0)
    graphPoints(fmatrix)
    
    strokeWeight(2)
    stroke(255,0,0)
    graphPoints(transpose(multmatrix(transformation_matrix,transpose(fmatrix))))
    
    
#用线段连接列表中的相邻的点
def graphPoints(matrix):
    beginShape()
    for pt in matrix:
        vertex(pt[0] * xscl,pt[1] * yscl)
    endShape()
​
def grid(xscl,yscl):
    strokeWeight(1)
    stroke(0,255,255)
    for i in range(xmin,xmax+1):
        line(i*xscl,ymin*yscl,i*xscl,ymax*yscl)
    for i in range(ymin,ymax+1):
        line(xmin*xscl,i*yscl,xmax*xscl,i*yscl)
    stroke(0)
    line(0,ymin*yscl,0,ymax*yscl)
    line(xmin*xscl,0,xmax*xscl,0)
    
#a为取行，b为取列
def multmatrix(a,b):
    m = len(a)
    n = len(b[0])
    
    newmatrix = []
    
    for i in range(m):
        row = []
        for j in range(n):
            sum1= 0
            for k in range(len(b)):
                sum1 += a[i][k] * b[k][j]
            row.append(sum1)
        newmatrix.append(row)
    return newmatrix
​
#转置函数
def transpose(a):
    newmatrix = []
    for i in range(len(a[0])):
        row = []
        for j in range(len(a)):
            row.append(a[j][i])
        newmatrix.append(row)
    return newmatrix