本文内容来自于学习麻省理工学院公开课：单变量微积分-分部积分-网易公开课
开发环境准备：CSDN
目录
一、多项式部分分式方法求积分
1、效果
2、步骤
(1) 长除法
(2) 分解因式 (factor)
(3) 建立等式（setup）
二、分部积分(integration by parts)
1、公式
2、例子
(1)
(2)
3、换算公式（Reduction formula）
4、
5、
6、求 绕 y 轴旋转， 形成的体的体积，
一、多项式部分分式方法求积分
多项式相除
1、效果(1) 总是可以用
(2) 处理起来有些慢
2、步骤 (1) 长除法 (2) 分解因式 (factor)这里
(3) 建立等式（setup）注意这里有12个未知数 , 要列12个等式来解决，而每个等式都会非常长，需要求助计算机 。
假设这个R(x) = x+1, 那么求 的积分
import numpy as np
from sympy import *
x = symbols('x')
expr = (x+1)/((x+2)**4*(x**2+2*x+3)*(x**2+4)**3)
print (float(integrate(expr, (x, 1, 5))))1.1996123519913639e-05
对已经分开的式子求积分
a.
b.
由公式：
设
由公式： (半角公式)
由公式：
(倍角公式)
(倍角公式)
由 ，
expr = 1/(x**2+4)**3
print(integrate(expr,x))(3*x**3 + 20*x)/(128*x**4 + 1024*x**2 + 2048) + 3*atan(x/2)/256
对于这个复杂度，老师是这么讲的，积分有时就是这么复杂，对谁来说都是这样，把它交给计算机去算就好，可以看到如果用计算机处理的话，就两行....
二、分部积分(integration by parts) 1、公式由
有
有公式：
2、例子 (1)设
则
(2)设
则
设
则
3、换算公式（Reduction formula）求
总结：
n-1, n-1->n-2, ....->0)" src="https://latex.codecogs.com/gif.latex?F_%7Bn%7D%28x%29%20%3Dxln%5E%7Bn%7D%28x%29-n%28F_%7Bn-1%7D%29%20%2C%20%28n-%3En-1%2C%20n-1-%3En-2%2C%20....-%3E0%29" />
4、设
则
5、设
则
当n=0
当n=1
当n=2
设
则
总结：
检查：
u= symbols('x')
expr1 = x**3*cos(x)
print ( 'f3 = ', integrate(expr1))
expr1 = x**4*cos(x)
print ( 'f4 = ',integrate(expr1))
x= symbols('x')
expr1 = x**3*sin(x) + 3*x**2*cos(x) - 6*x*sin(x) - 6*cos(x)
print('d(f3) = ', diff(expr1))f3 = x**3*sin(x) + 3*x**2*cos(x) - 6*x*sin(x) - 6*cos(x)
f4 = x**4*sin(x) + 4*x**3*cos(x) - 12*x**2*sin(x) - 24*x*cos(x) + 24*sin(x)
d(f3) = x**3*cos(x)
这里用程序计算了 ,如上
检查结果：符合
6、求 绕 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(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 = np.e**(x)
xx,yy,zz = MakeRotateByAxisS(0,1,50,expr,0.0,2*np.pi,50, EnumAxis.YAxis)
points = GridToArray(xx,yy,zz)
ax.plot_surface(xx,yy,zz)
ax.view_init(elev=40, # 仰角
azim=110 # 方位角
)
plt.show()
圆盘法：
Bullseye：第三单元 用python学习微积分（二十）壳层法、圆盘法求体积 (上）
由公式：
壳层法:
Bullseye：第三单元 用python学习微积分（二十）壳层法、圆盘法求体积 (下）
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）