第四单元用python学习微积分（二十八）参数方程、弧长和表面积

本文内容来自于学习麻省理工学院公开课：单变量微积分-分部积分-网易公开课Bullseye：第一单元用python学习微积分（一）安装开发环境Anaconda 和导数（上）- 1/x的导数

Bullseye

1858人浏览 · 2022-05-02 17:26:40

Bullseye · 2022-05-02 17:26:40 发布

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

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

5、例2 -- 半圆的弧长 ( 弧长和sin的定义推导 )

2、球面的表面积(Surface Area of A Sphere)

三、参数方程（Parametric Curves)

一、弧长

1、前提

曲线 S 经过点 $(s_0,s_1...,s_n)$ , 这些点均匀分布，并分别投影在x轴上 $(a,x_1,...b)$

当计算曲线 S 的长度时，我们可以通过用直线连接这些点 $(s_0,s_1...,s_n)$ ，并且累加这些线段的长度以求得近似解。

2、对一小段曲线进行长度计算：

每一小段弧长的长度

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

而微积分的思想是，当在极限的情况下，也就是当曲线被分为无限多的份数时， $\Delta x, \Delta y$ 无限小, 有

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

（注意：这里老师提到上式可以简写作 $ds^2 = dx^2 + dy^2$ ）

(1) 公式1： $ds = \sqrt{dx^2+dy^2}$

(2) 公式2： $ds = \sqrt{1+(\frac{dy}{dx})^2} \times dx$

3、总弧长公式 $(s_n-s_0)$

$s = \int_a^b \sqrt{1+(\frac{dy}{dx})^2} \times dx (=\int ds)$

$= \int_a^b \sqrt{1+f'(x)^2} \times dx....... (y = f(x))$

4、例1-直线的长度

$y = mx$

$y' = m$

$ds = \sqrt{1+(y')^2}dx = \sqrt{1+m^2}dx$

曲线在x轴上的投影 $0 \leq x \leq10$

$s = \int_0^{10} \sqrt{1+m^2}dx = \sqrt{1+m^2} x|_0^{10} = 10\sqrt{1+m^2}$

import numpy as np 
from sympy import *
import matplotlib.pyplot as plt 
figure, ax= plt.subplots( 1 ) 
ax.set_aspect( 1 ) 
def DrawAngleSign(radiansFrom, radiansTo, radius, axes, color):
    angleInner = np.linspace(  radiansTo , radiansFrom, 150 ) 
    xArcInner = radius /10 * np.cos( angleInner ) 
    yArcInner = radius /10 * np.sin( angleInner ) 
    axes.plot( xArcInner, yArcInner,color=color) 

def DrawXY(xFrom,xTo,steps,expr,color,label,plt):
    yarr = []
    xarr = np.linspace(xFrom ,xTo, steps) 
    for xval in xarr:
        yval = expr.subs(x,xval)
        yarr.append(yval)
    y_nparr = np.array(yarr) 
    plt.plot(xarr, y_nparr, c=color, label=label)  
    
m = 0.5
x = symbols('x')
expr = m*x
DrawXY (0,10, 50, expr, color='c',label='y = mx',plt=plt)
plt.plot([10,10], [0,5], c='g')  
plt.plot([0,10], [0,0], c='g')  
plt.legend(loc='upper left')
plt.show()

y=mx

这条蓝色的线长： $\sqrt{x^2+y^2} = \sqrt{10^2 + m10^2} = 10\sqrt{1+m^2}$

5、例2 -- 半圆的弧长 ( 弧长和sin的定义推导 )

$y = \sqrt{1-x^2}$

$y' = \frac{1}{2}(1-x^2)^{-\frac{1}{2}} (-2x) = -x(1-x^2)^{-\frac{1}{2}}$

$ds = \sqrt{1+(y')^2}dx = \sqrt{1+(-x(1-x^2)^{-\frac{1}{2}})^2}dx = \sqrt{1+\frac{x^2}{1-x^2}} = \sqrt{\frac{1}{1-x^2}}$

$\alpha = \int_0^a \sqrt{\frac{1}{1-x^2}}dx = arcsin(x)|_0^a = arcsin(a)$

$\alpha = arcsin(a)$

$sin(\alpha) = a$

import numpy as np 
import matplotlib.pyplot as plt 
 
angle = np.linspace( 0 , 2 * np.pi , 150 ) 
 
radius = 1
 
x = radius * np.cos( angle ) 
y = radius * np.sin( angle ) 
 
figure, axes = plt.subplots( 1 ) 
 
axes.plot( x, y, label='circle radius:'+format(radius),color='b' ) 


axes.set_aspect( 1 ) 
x = symbols('x')
expr = (1-x**2)**0.5
a = 0.5

DrawXY(0,0.5,50,expr,color='c',label='y = (1-x**2)**0.5',plt=plt)

plt.plot([0,0], [0,1], c='g')  
plt.plot([0.5,0.5], [0,expr.subs(x, 0.5)], c='red',label='x=a')  
plt.plot([0,0.5], [0,expr.subs(x, 0.5)],c='g')  
plt.plot([0,1], [0,0],c='b')  
plt.text(-0.2, -0.1, 'sin(alpha) = a ', fontsize=12)
plt.text(0.1, 1.0, 'alpha', fontsize=12)
plt.plot([0.5],[0],lw=0, marker='o', color='red',fillstyle='none')
plt.text(0.03, 0.3, 'alpha', fontsize=9)

alhpa = asin(a)
DrawAngleSign(0.5*np.pi-float(alhpa), 0.5*np.pi, 1,axes, 'b')


plt.title( 'Circle' ) 
plt.legend(loc='lower left')
plt.show()

$sin(alpha) = a$

6、例3 -- 抛物线的长度

$y = x^2$

$y' = 2x$

$ds = \sqrt{1+(y')^2}dx = \sqrt{1+4x^2}dx$

$0\leq x\leq a$

$\int_0^a \sqrt{1+4x^2}dx$

由三角替换公式：

$\sqrt{a^2 +x^2} -- x = a \times tan(\theta) -- a \times sec(\theta) -- dx = asec^2(\theta)d\theta$

设 $x =\frac{1}{2} tan(u)$

则 $dx = \frac{1}{2}sec^2(u)du$

$\int\sqrt{1+4x^2}dx = \frac{1}{2}\int sec^3(u)du$

这里还是计算机吧....

sympy无法计算结果，因为它没有找到这个积分的原函数

只能使用数值积分得到结果，如计算 $\int_0^5 \sqrt{1+4x^2}dx$

from scipy import integrate
result=integrate.quad(lambda x: (1+4*x**2)**0.5, 0,5)
print(result)

(25.87424479037674, 6.135691242135164e-08)

结果是25.87424479037674...

二、表面积

1、抛物线旋转成面的表面积

有曲线 $y=x^2$ ，围绕 x 轴转一周

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(xFrom,xTo,steps,expr,thetaFrom,thetaTo,thetaSteps, rotatedBy):
    stepsArr = np.linspace(thetaFrom ,thetaTo, thetaSteps)
    xarr = []
    yarr = []
    zarr = []
    i = 0
    for theta in stepsArr:
        x,y,z = MakeRotateByAxis(xFrom,xTo,steps,expr,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(xFrom,xTo,steps,expr,theta,rotatedBy):
    yarr = []
    xarr = np.linspace(xFrom ,xTo, steps) 
    xyz = []
    xx = []
    yy = []
    zz = []
    for xval in xarr:
        yval = expr.subs(x,xval)
        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(xFrom,xTo,steps,expr,theta,color,ax,rotatedBy):
    xx,yy,zz = MakeRotateByAxis(xFrom,xTo,steps,expr,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)

x = symbols('x')
expr = x**2
xx,yy,zz = MakeRotateByAxisS(0,1,50,expr,0.0,2*np.pi,50, EnumAxis.XAxis)
points = GridToArray(xx,yy,zz)
ax.plot_surface(xx,yy,zz)
ax.view_init(elev=40,    # 仰角
             azim=110    # 方位角
            )
plt.show()

抛物线绕 x 轴旋转

(1) 表面积计算公式，类壳层法

求表面积

考虑上图中的 ds 很段，则 dA 相当于求一个圆柱的表面积（不计算顶、底）

这里老师介绍道：一般用 s 代表弧长，而 S 代表表面积

这里 $dA = (2\pi y)(ds)$ 也可以写成 $dS = (2\pi y)(ds)$

$y = x^2$

由抛物线的长度的公式

注意这里计算ds是在x，y平面，所以适用弧长公式

$ds = \sqrt{1+(y')^2}dx =\sqrt{1+4x^2}dx$

$Surface Area = \int_0^a 2\pi yds = \int_0^a 2\pi x^2 \sqrt{1+4x^2}dx$

from scipy import integrate
result=integrate.quad(lambda x: (2*np.pi*x**2*(1+4*x**2)**0.5)**0.5, 0,5)
print(result)

(80.00795040596235, 1.1729232113152397e-08)

2、球面的表面积(Surface Area of A Sphere)

半径: a

$y = \sqrt{a^2 - x^2}$

$y' = \frac{-2x}{2\sqrt{a^2-x^2}} = \frac{-x}{\sqrt{a^2-x^2}}$

注意这里计算ds是在x，y平面，所以适用弧长公式

由公式： $ds = \sqrt{1+(y')^2}dx = \sqrt{\frac{a^2}{a^2-x^2}}dx= a\sqrt{\frac{1}{a^2-x^2}}dx$

在 x, y 平面画一个半圆，然后绕 x 轴旋转360度，形成球体

fig = plt.figure(figsize=(12, 12),
                facecolor='lightyellow'
                )
# draw sphere
ax = fig.gca(fc='whitesmoke',
              projection='3d' 
           )
DrawAxis(ax, 3,3,3)

x = symbols('x')
expr = (1 - x**2)**0.5
xx,yy,zz = MakeRotateByAxisS(-1,1,50,expr,0.0,2*np.pi,50, EnumAxis.XAxis)
points = GridToArray(xx,yy,zz)
ax.plot_surface(xx,yy,zz)
ax.view_init(elev=40,    # 仰角
             azim=110    # 方位角
            )
plt.show()

球

要得到球的表面积 $a\leq x\leq b$ ,设 a=-0.1, b=0.1, 把上面代码改这一行并执行

xx,yy,zz = MakeRotateByAxisS(-0.1,0.1,50,expr,0.0,2*np.pi,50, EnumAxis.XAxis)

-0.1<x<0.1

可以看出还是可以用上例的方法计算，也就是把这部分球看成一个个高度很小圆柱拼接成，然后计算。

$Surface Area = \int_{x_1}^{x_2} 2\pi yds = \int_{x_1}^{x_2} 2\pi a\sqrt{a^2-x^2} \sqrt{\frac{1}{a^2-x^2}}dx= \int_{x_1}^{x_2} 2\pi adx$

$=2\pi a(x_2-x_1)$

验证：半径为a

半球的表面积： $0\leq x \leq a , A = 2\pi a(a-0) = 2\pi a^2$

球的表面积： $-a\leq x \leq a , A = 2\pi a(a-(-a)) = 4\pi a^2$

符合球的表面积公式

三、参数方程（Parametric Curves)

1、前提

t 是参数 ( parameter )

$x = x(t)$

$y = y(t)$

举例：t表示时间，某物体沿轨道 ( trajectory )行进，当 $t = t_1$ 时，物体的位置为 ( $x(t_1), y(t_1)$ )

物体沿轨道行进

2、例1 -- ⚪

(1) 形状

$x = acos(t)$

$y = asin(t)$

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

从这里可以看出这个由 { x, y, t } 共同组成的方程组描绘的是一个半径为a的⚪

(2) 方向

$t = 0, (x, y) = (acos(0), asin(0)) = (0,a)$

$t =\frac{\pi}{2}, (x, y) = (acos(\frac{\pi}{2}), asin(\frac{\pi}{2})) = (a,0)$

...

由以上可以看出，以当 t 增长时，函数的方向是逆时针的（counter clockwise）

(3) 弧长

当长度很短时，可以把线段长度看作是弧长

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

验证：

周长： $\int_0^{2\pi}adt = at|_0^{2\pi} = 2\pi a$

腾讯云开发者社区

腾讯云面向开发者汇聚海量精品云计算使用和开发经验，营造开放的云计算技术生态圈。

更多推荐

动物识别系统 python实现+UI

腾讯云开发者社区

Linux下安装MySQL8.0(超详细)

腾讯云开发者社区

安卓手机如何连接苹果电脑

腾讯云开发者社区

所有评论(0)

查看更多评论

Bullseye

@bullseye

已为社区贡献2条内容

第四单元 用python学习微积分（二十八）参数方程、弧长和表面积

Bullseye

一、弧长

1、前提

2、对一小段曲线进行长度计算：

(1) 公式1：

(2) 公式2：

3、总弧长公式

4、例1-直线的长度

5、例2 -- 半圆的弧长 ( 弧长和sin的定义推导 )

6、例3 -- 抛物线的长度

二、表面积

1、抛物线旋转成面的表面积

(1) 表面积计算公式，类壳层法

2、球面的表面积(Surface Area of A Sphere)

三、参数方程（Parametric Curves)

1、前提

2、例1 -- ⚪

(1) 形状

(2) 方向

(3) 弧长

所有评论(0)

Bullseye

第四单元用python学习微积分（二十八）参数方程、弧长和表面积

(1) 公式1： $ds = \sqrt{dx^2+dy^2}$

(2) 公式2： $ds = \sqrt{1+(\frac{dy}{dx})^2} \times dx$

3、总弧长公式 $(s_n-s_0)$