第四单元用python学习微积分（二十九）参数曲线(二)和极坐标和极坐标下的面积

本文链接：https://blog.csdn.net/bullseye/article/details/125036813

本文内容来自于学习麻省理工学院公开课：单变量微积分-参数方程、弧长和表面积-网易公开课

Bullseye：第一单元用python学习微积分（一）安装开发环境Anaconda 和导数（上）- 1/x的导数

Bullseye：第一单元用python学习微积分（二）VSCode 、PYGame 和导数（上）- 瞬时速度

(1) 之前举例的曲线的参数方程共同的前提条件：

3、椭圆 x = 2 sin(t); y = cos(t)

(1) 参数方程

4、椭球的表面积

二、极坐标（polar coordinate）

3、例1 点 ( 1, -1 ) 用极坐标表示

(1)

(2)

(3) 这里可以看成是原点在方向的反向延申的距离

一、参数曲线（续）

1、⚪

$x = acos(t)$

$y = asin(t)$

$x^2+ y^2 = a^2cos^2(t) +a^2sin^2(t) = a^2$

以当 t 增长时，函数的方向是逆时针的（counter clockwise）

(1) 求弧长：

$ds^2 = dx^2+dy^2$

$ds = \sqrt{(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2} dt$

$\frac{dx}{dt} = -asin(t)$

$\frac{dy}{dt} = acos(t)$

(2) 速度

$ds = \sqrt{a^2sin^2(t) + a^2cos^2(t)}dt = adt$

$\frac{ds}{dt} = a$ ...... (可以考虑a是速度，围绕圆心，以a为速度环绕)

(3) 修改速度

$x = acos(kt)$

$y = asin(kt)$

import numpy as np 
from sympy import *
a,k,x = symbols('a,k,x')
expr = a*sin(k*x)
print('a*sin(kx)=', diff(expr,x))
expr = a*cos(k*x)
print('a*cos(kx)=', diff(expr,x))

a*sin(kx)= a*k*cos(k*x) a*cos(kx)= -a*k*sin(k*x)

$ds = \sqrt{(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2} dt = \sqrt{a^2k^2cos^2(kx) +a^2 k^2sin^2(kx)} dt=akdt$

$\frac{ds}{dt} = ak$ ...... ( 这时ak是新的速率 )

( a>0, k>0 )

2、注意：

(1) 之前举例的曲线的参数方程共同的前提条件：

$\Delta s \approx \Delta x^2 + \Delta y^2$

$(\frac{ds}{dt})^2 = (\frac{dx}{dt})^2+(\frac{dy}{dt})^2$

$\frac{ds}{dt} = \sqrt{(\frac{dx}{dt})^2 + (\frac{dy}{dt})^2}$

这里老师举了个反例： $\frac{dx^2}{dt^2}$ ！！这种写法是不容许的！！

3、椭圆 x = 2 sin(t); y = cos(t)

(1) 参数方程

$\frac{1}{4} x^2 + y^2 = sin^2t + cos^2t = 1$ （方向顺时针）

$\frac{ds}{dt} = \sqrt{(\frac{dx}{dt})^2+(\frac{dy}{dt})^2} = \sqrt{(2cos(t))^2+sin^2(t)}$

弧长= $\int_{0}^{2\pi}\sqrt{4cos^2(t)+sin^2(t)}dt$ (老师语：这个不是初等积分，到此处就是答案)

对这个弧长做个数值积分：

from scipy import integrate
result=integrate.quad(lambda x: (4*cos(x)**2+ sin(x)**2)**0.5, 0,2*pi)
print(result)
(9.688448220547677, 3.467269449700442e-10)

这是个椭圆，注意 x,y 和 t 的关系。

我在第二讲时写了个例子用到了pygame，这次还是会用到它，进行了简单的重构，添加了几个新类，增加了画线和画圆的功能，改掉了几个bug。代码打包在百度网盘请输入提取码提取码：ywmc

关于代码，说点题外话，简单说说重构。这次的需求和上次相比，多了画圆形和画线的功能，可以看出它们都是图形，最好有一个共同的基类，类似的操作封装在一个函数里就好，这里是考验思维得地方，多想想怎么做才是最节省，最清晰的。没花太多时间的修改，还是不太好，但是很高兴的是原本很简陋的代码经过添加了几个功能变得好了一点点。在开发中常会遇到类似的问题，开发个功能相似的模块怎么弄？如果去拷贝、黏贴、改一改，每个类都有自己完整的一套api，那类似的代码会被贴得到处都是，值得注意的是随之拷贝的还有旧代码中原生的错误，换句话说就是错误会被拷贝的到处都是。想象一下，这之后有人在某个功能上发现了那个错误提交修改，而负责修改的人或许只看到了这一处，那他就对症改改好，修改中也有可能带来其他的错误，而后又有人去拷贝他的，周而复始，传说中的屎山就是这么来的。

所以追求代码的质量，减少冗余，其实就是在追求系统的健壮性和减少你的工作量。不需要求一步到位，只要每次改代码或添加功能时，想着让这段代码变好一点，抱着这种想法，你也许就会收获更轻松的工作。

主程序代码：

import random
import pygame
import math
import sys
import time

from pygame.constants import K_BACKSPACE, K_END, K_RETURN
from Coordinate import *
from screen import *
from XRectangle import XRectangle
from graphic import *
from TextBox import *
from numpy import *

def CalculateDistance1(seconds, speed):
    return speed*cos(seconds), speed*sin(seconds)

def CalculateDistance2(seconds, speedX, speedY):
    return speedX*sin(seconds), speedY*cos(seconds)

def IsPtInList(pt, ptList):
    for i in ptList:
        if(i[0] == pt[0] and i[1] == pt[1]):
            return True
    return False
    
if __name__ == "__main__":
    x0 = 300
    y0 = 300
    height = 2
    r = 10
    speedY = speedX = 100
    timesOfSpeed = 1
    scr = Screen()
    clock = pygame.time.Clock()

    pygame.init()
    coordinate = Coordinate(scr, 1, False)
    pointOnCircle = Circle(x0, y0, r*2, r*2, coordinate, scr, 1)
    
    centerPoint = Circle(x0, y0, r*2, r*2, coordinate, scr, 1)
    
    enterInfo = TextBox(100, 590,24,coordinate,scr)
    tabInfo = TextBox(100, 560,24,coordinate,scr)
    speedXInfo = TextBox(100, 530,24,coordinate,scr)
    timeInfo = TextBox(100, 500,24,coordinate,scr)
    
    clock = pygame.time.Clock()
    timeBegin = 0
    seconds = 0
    running = True
    xx = 0
    tmpDistY = 0
    tmpDistY1 = 0
    tmpSpeed = '0'
    tmpSeconds = 0
    dist = 0
    timeBegin = time.time()
    lsPoints = []
    while running:
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                running = False
                break
            elif event.type == pygame.KEYDOWN:
                key = event.key
                if key == K_RETURN:
                    lsPoints.clear()
                    if(timesOfSpeed!=3):
                        timesOfSpeed = 3
                    else:
                        timesOfSpeed = 1
                elif key == K_BACKSPACE:
                    lsPoints.clear()
                    timeBegin = time.time()
                elif key == K_END:
                    lsPoints.clear()
                    if(speedY == speedX):
                        speedX = 2 * speedY
                    else:
                        speedX = speedY
        scr.fill()
        centerPoint.draw()
        if timeBegin!=0:
            seconds = time.time()
            seconds = seconds - timeBegin
            if(speedY == speedX):
                distX,distY = CalculateDistance1(seconds*timesOfSpeed, speedX)
            else:
                distX,distY = CalculateDistance2(seconds*timesOfSpeed, speedX,speedY)
            
            pointOnCircle.center()
            pointOnCircle.moveTo(x0+distX, y0+distY)
            pointOnCircle.draw()
            centerPoint.draw()
            center0 = (x0+distX , y0+distY)
            center0 = pointOnCircle.GetPointData(center0)
            if(not IsPtInList(center0, lsPoints)):
                lsPoints.append(center0)
            pointOnCircle.drawLine(lsPoints)
        else:
            centerPoint.draw()
            pointOnCircle.drawLine(lsPoints)
            
        timeInfo.draw('timer:' +str(seconds))
        tabInfo.draw('press backspace to restart')
        if(timesOfSpeed == 1):
            enterInfo.draw('press enter to make it 3 times faster')
        else:
            enterInfo.draw('press enter to make it back to normal')
        speedXInfo.draw("press end to change formula")
        pygame.display.flip()
        clock.tick(100)

4、椭球的表面积

上面椭圆沿y轴旋转一周，形成的椭球

from sympy import *
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from enum import Enum
from itertools import product, combinations

from matplotlib import cm
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(12, 12),
                facecolor='lightyellow'
                )
# draw sphere
ax = fig.gca(fc='whitesmoke',
              projection='3d' 
           )
 
class EnumAxis(Enum):
    XAxis = 1
    YAxis = 2
    ZAxis = 3
    
def fromXYZM(xyzM):
    ret = []
    
    m=range(0,xyzM.shape[0])
    
    for b in zip(m,[0]):
        ret.append((xyzM[b]))
    for b in zip(m,[1]):
        ret.append((xyzM[b]))
    for b in zip(m,[2]):
        ret.append((xyzM[b]))
    return ret

def plot_surface(x, y, z, dx, dy, dz, ax):
    xx = np.linspace(x, x+dx, 2)
    yy = np.linspace(y, y+dy, 2)
    zz = np.linspace(z, z+dz, 2)
    
    if(dz == 0):
        xx, yy = np.meshgrid(xx, yy)
        ax.plot_surface(xx, yy, z, alpha=0.25)
    
    if(dx == 0):
        yy, zz = np.meshgrid(yy, zz)
        ax.plot_surface(x, yy, zz, alpha=0.25)
    
    if(dy == 0):
        xx, zz = np.meshgrid(xx, zz)
        ax.plot_surface(xx, y, zz, alpha=0.25)
    

def plot_opaque_cube(x, y, z, dx, dy, dz, ax):

    xx = np.linspace(x, x+dx, 2)
    yy = np.linspace(y, y+dy, 2)
    zz = np.linspace(z, z+dz, 2)

    xx, yy = np.meshgrid(xx, yy)
    ax.plot_surface(xx, yy, z)
    ax.plot_surface(xx, yy, z+dz)

    yy, zz = np.meshgrid(yy, zz)
    ax.plot_surface(x, yy, zz)
    ax.plot_surface(x+dx, yy, zz)

    xx, zz = np.meshgrid(xx, zz)
    ax.plot_surface(xx, y, zz)
    ax.plot_surface(xx, y+dy, zz)
    # ax.set_xlim3d(-dx, dx*2, 20)
    # ax.set_xlim3d(-dx, dx*2, 20)
    # ax.set_xlim3d(-dx, dx*2, 20)
    
def DrawAxis(ax, xMax, yMax, zMax):
    # 设置坐标轴标题和刻度
    ax.set(xlabel='X',
           ylabel='Y',
           zlabel='Z',
           xticks=np.arange(-xMax, xMax, 1),
           yticks=np.arange(-yMax, yMax, 1),
           zticks=np.arange(-zMax, zMax, 1)
          )

    ax.plot3D(xs=[xMax,-xMax, 0,0, 0, 0, 0,0, ],    # x 轴坐标
              ys=[0, 0,0, yMax,-yMax, 0, 0,0, ],    # y 轴坐标
              zs=[0, 0,0, 0,0, 0, zMax,-zMax ],    # z 轴坐标
              zdir='z',    # 
              c='k',    # color
              marker='o',    # 标记点符号
              mfc='r',    # marker facecolor
              mec='g',    # marker edgecolor
              ms=10,    # size
            )

def Rx(x,y,z,theta):
    A = [x,y,z] * np.matrix([[ 1, 0           , 0           ],
                   [ 0, cos(theta),-sin(theta)],
                   [ 0, sin(theta), cos(theta)]])
    return fromXYZM(A)
 
def Ry(x,y,z,theta):
    A = [x,y,z] * np.matrix([[ cos(theta), 0, sin(theta)],
                   [ 0           , 1, 0           ],
                   [-sin(theta), 0, cos(theta)]])
    return fromXYZM(A)
 
def Rz(x,y,z,theta):
    A = [x,y,z] * np.matrix([[ cos(theta), -sin(theta), 0 ],
                   [ sin(theta), cos(theta) , 0 ],
                   [ 0           , 0            , 1 ]])
    return fromXYZM(A)

def Equal(a, b, tol):
    if(abs(a - b) < tol):
        return true
    return false

def GetValuesByIdxMGrid(grid, idx):
    length = int(len(idx)/2)
    values = []
    for n in range(length):
        i = idx[2*n]
        j = idx[2*n+1]
        values.append(grid[i][j])
    return values
        
def FindIndexInMGrid(grid, value, tol):
    idx = []
    values = []
    for i in range(len(grid)):
        for j in range(len(grid[i])):
            if(Equal(grid[i][j], value, tol)):
                idx.append(i)
                idx.append(j)
                values.append(value)
                
    return idx,values
            
def MakeSliceByAxis(xarr, yarr, zarr, axis, value, tol):
    idx = []
    xret = []
    yret = []
    zret = []
    if axis == EnumAxis.XAxis:
        idx,xret = FindIndexInMGrid(xarr, value,tol)
        yret = GetValuesByIdxMGrid(yarr, idx)
        zret = GetValuesByIdxMGrid(zarr, idx)
    else:
        if axis == EnumAxis.YAxis:
            idx,yret = FindIndexInMGrid(yarr, value,tol)
            xret = GetValuesByIdxMGrid(xarr, idx)
            zret = GetValuesByIdxMGrid(zarr, idx)
        else:
            idx,zret = FindIndexInMGrid(zarr, value,tol)
            xret = GetValuesByIdxMGrid(xarr, idx)
            yret = GetValuesByIdxMGrid(yarr, idx)
            
            
    return np.array(xret),np.array(yret),np.array(zret)
    
    
def MakeRotateByAxisS(tFrom,tTo,steps,exprX, exprY,thetaFrom,thetaTo,thetaSteps, rotatedBy):
    stepsArr = np.linspace(thetaFrom ,thetaTo, thetaSteps)
    xarr = []
    yarr = []
    zarr = []
    i = 0
    for theta in stepsArr:
        x,y,z = MakeRotateByAxis(tFrom,tTo,steps,exprX,exprY, theta,rotatedBy)
        xarr.insert(i, x)
        yarr.insert(i, y)
        zarr.insert(i, z)
        i = i + 1
    
    return np.array(xarr), np.array(yarr), np.array(zarr)

def MakeRotateByAxis(tFrom,tTo,steps,exprX,exprY, theta,rotatedBy):
    yarr = []
    xarr = []
    tarr = np.linspace(tFrom ,tTo, steps) 
    xyz = []
    xx = []
    yy = []
    zz = []
    for tval in tarr:
        yval = exprY.subs(t,tval)
        xval = exprX.subs(t,tval)
        if rotatedBy == EnumAxis.XAxis:
            xyz =  Rx(xval,yval,0,theta)
        elif rotatedBy == EnumAxis.YAxis:
            xyz =  Ry(xval,yval,0,theta)
        else:
            xyz =  Rz(xval,yval,0,theta)
            
        xx.append(float(xyz[0])) #using float() to prevent isnan issue
        yy.append(float(xyz[1]))
        zz.append(float(xyz[2]))
        #if(np.isnan(float(xyz[2]))):
            #zz.append(float(0.0))
        #else:
            #zz.append(xyz[2])
    return xx,yy,zz
    
def RotateByAxis(tFrom,tTo,steps,exprX,exprY, theta,color,ax,rotatedBy):
    xx,yy,zz = MakeRotateByAxis(tFrom,tTo,steps,exprX,exprY, theta,rotatedBy)
    ax.plot(xx, yy, zz, color)
    
        
def GridToArray(xx,yy,zz):
    ret = []
    length = len(xx)
    n = 0
    for i in range(length):
        length1 = len(xx[i])
        for j in range(length1):
            coordinate = []
            coordinate.append(xx[i][j])
            coordinate.append(yy[i][j])
            coordinate.append(zz[i][j])
            ret.insert(n, coordinate)
            n += 1
    return np.array(ret) 
                   


DrawAxis(ax, 3,3,3)

t = symbols('t')
exprX = 2*sin(t)
exprY = cos(t)

xx,yy,zz = MakeRotateByAxisS(0,pi,50,exprX,exprY, 0.0,2*np.pi,50, EnumAxis.YAxis)
points = GridToArray(xx,yy,zz)
ax.plot_surface(xx,yy,zz, color='green', alpha=0.45)
ax.view_init(elev=50,    # 仰角
             azim=40    # 方位角
            )
plt.show()

求这个体的表面积，还是得从x，y坐标系考虑。因为这是一个在x，y平面的⚪，绕y轴旋转生成的体。

任意一点的 $x= 2sin(t)$ ；

所对应的绕y轴的⚪的周长为 $2\pi x= 2\pi sin(t)$

这个⚪周长乘以一个极小的弧长（表示高差）可以计算出这个这个⚪为底高差为高的环形的表面积，累加（积分）即可得到这个椭球的表面积

由上小节结论可知弧长= $\int_{0}^{2\pi}\sqrt{4cos^2(t)+sin^2(t)}dt$

表面积为 $da = \int_{0}^{\pi} 2\pi sin(t)\sqrt{4cos^2(t)+sin^2(t)}dt$

二、极坐标（polar coordinate）

1、定义

这是另一种坐标系，对平面上的一点的描述为一对参数 $( r, \theta )$

“r 是这个点到原点（坐标系两个轴的交点）的距离

$\theta$ 是这个点到原点的线段与水平轴的夹角”

(1) 准确的定义：

$x = r cos(\theta)$

$y =rsin(\theta)$

(2) 推导出的公式：

$r = \sqrt{x^2+y^2}$

$\theta = tan^{-1}\frac{y}{x}$ ...... $(arctan(\frac{y}{x}))$

这两个式子有模糊的地方，有以下这些情况，需要看图决定

$r = \pm\sqrt{x^2+y^2}$

$\theta = tan^{-1}\frac{-y}{-x}$

2、坐标系介绍

(1) 直角坐标系

plt.plot([-2,-2,-2], [-2,0,2], c='r', label='x=-2')  
plt.plot([0,0,0], [-2,0,2], c='r', label='x=0')  
plt.plot([2,2,2], [-2,0,2], c='r', label='x=2')  

plt.plot([-2,0,2], [-2,-2,-2], c='c', label='y=-2')  
plt.plot([-2,0,2], [0,0,0], c='c', label='y=0')  
plt.plot([-2,0,2], [2,2,2], c='c', label='y=2')  
plt.legend(loc='upper left')
plt.show()

(2) 极坐标系

import numpy as np 
from sympy import *
import matplotlib.pyplot as plt 
figure, ax= plt.subplots( 1 ) 
ax.set_aspect( 1 ) 


def DrawXY(tFrom,tTo,steps,exprX,exprY, color,label,plt):
    xarr = []
    yarr = []
    tarr = np.linspace(tFrom ,tTo, steps) 
    for tval in tarr:
        xval = exprX.subs(t,tval)
        xarr.append(xval)
        yval = exprY.subs(t,tval)
        yarr.append(yval)
    y_nparr = np.array(yarr) 
    x_nparr = np.array(xarr) 
    plt.plot(x_nparr, y_nparr, c=color, label=label)  

r = 1    
t = symbols('t')
exprX = r*cos(t)
exprY = r*sin(t)
DrawXY( 0,2*np.pi,50,exprX,exprY,color='c', label='r=1',plt = plt)
r = 2
exprX = r*cos(t)
exprY = r*sin(t)
DrawXY( 0,2*np.pi,50,exprX,exprY,color='g', label='r=2',plt = plt)

a = 0
exprX = t*cos(a)
exprY = t*sin(a)
DrawXY( -2,2,50,exprX,exprY,color='blue', label='t=0',plt = plt)

a = pi/4
exprX = t*cos(a)
exprY = t*sin(a)
DrawXY( 0,2,50,exprX,exprY,color='red', label='t=pi/4',plt = plt)

a = pi/4
exprX = t*cos(a)
exprY = t*sin(a)
DrawXY( -2,0,50,exprX,exprY,color='orange', label='t=-pi/4',plt = plt)

plt.legend(loc='upper left')
plt.show()

3、例1 点 ( 1, -1 ) 用极坐标表示

$r = \pm\sqrt{x^2+y^2}$

$r = \pm\sqrt{2}$

(1) $(r,\theta) = (\sqrt{2},\frac{7\pi}{4})$

(2) $(r,\theta) = (\sqrt{2},\frac{-\pi}{4})$

(3) $(r,\theta) = (-\sqrt{2},\frac{3\pi}{4})$ 这里可以看成是原点在 $\frac{3\pi}{4}$ 方向的反向延申 $\sqrt{2}$ 的距离

(4) $(r,\theta) = (-\sqrt{2},\frac{-5\pi}{4})$

4、例2 r = a 是圆形

5、例3 theta = c 是射线

这里注意，一般这里暗含了 $0<r\leq \infty$ ，当容许 $-\infty < r < \infty$ ,方程会给出一条直线

6、注意：使用极坐标常见的条件设置

$0\leq r <\infty$

$-\pi < \theta \leq \pi or 0 \leq \theta <2\pi$

7、例4 'y = 1'

$y = rsin(\theta)$

$r = \frac{1}{sin(\theta)} ( 0<\theta <\pi )$

t = symbols('t')
exprX = 1/sin(t)*cos(t)
exprY = 1/sin(t)*sin(t)
DrawXY( 0.1,np.pi-0.1,50,exprX,exprY,color='c', label='r=1/sin(t)',plt = plt)

8、例5 偏离原点的⚪

$(x-a)^2 +y^2 =a^2$

$x^2-2ax+a^2 + y^2 = a^2$

由： $x^2 + y^2 = r^2$

$r^2 -2ax = 0$

$r^2-2arcos(\theta) = 0$

$r = 2acos(\theta)$ or r = 0 .... ( $-\frac{\pi}{2} < \theta < \frac{\pi}{2}$ )

这里按老师的说法，和例四类似这里这个r是一个扫描的半径，在( $-\frac{\pi}{2} < \theta < \frac{\pi}{2}$ )区间把这个偏移出圆心的小圆的每个点扫出来

由：

$x = r cos(\theta)$

$y =rsin(\theta)$

figure, ax= plt.subplots( 1 ) 
ax.set_aspect( 1 ) 
a = 2
t = symbols('t')
r = 2*a*cos(t)
exprX = r*cos(t)
exprY = r*sin(t)
DrawXY( -np.pi/2,np.pi/2,50,exprX,exprY,color='c', label='r = 2*a*cos(t)',plt = plt)