栏目分类:
子分类:
返回
名师互学网用户登录
快速导航关闭
当前搜索
当前分类
子分类
实用工具
热门搜索
名师互学网 > IT > 软件开发 > 后端开发 > Python

第四单元 用python学习微积分(二十七)积分-部分分式-分部积分

Python 更新时间: 发布时间: IT归档 最新发布 模块sitemap 名妆网 法律咨询 聚返吧 英语巴士网 伯小乐 网商动力

第四单元 用python学习微积分(二十七)积分-部分分式-分部积分

本文内容来自于学习麻省理工学院公开课:单变量微积分-分部积分-网易公开课

开发环境准备: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学习微积分(二十)壳层法、圆盘法求体积 (下)

转载请注明:文章转载自 www.mshxw.com
本文地址:https://www.mshxw.com/it/850145.html
我们一直用心在做
关于我们 文章归档 网站地图 联系我们

版权所有 (c)2021-2022 MSHXW.COM

ICP备案号:晋ICP备2021003244-6号