python之微积分

ゝ一世哀愁。 2022-06-01 12:40 225阅读 0赞

**一、微分**

微分的话其实比较简单，手工算其实可以的，但是遇到特别复杂的函数，还是计算机算比较省时间，用到了sympy，举一个例子吧

import sympy as sy
    x=sy.symbols('x')#约定变量x
    y=x**3+10+sy.sin(x)#这个sin是sy的sin
    dy_dx=sy.diff(y,x)#常微分，写成dy_dx=sy.diff(y)也可以
    t=sy.symbols('t')
    z=x*sy.ln(t)+t**4
    dz_dt=sy.diff(z,t)#偏微分
    print(dz_dt)
    print(dy_dx)
    print(dy_dx.subs({x:5}))#把x=5代入
    print(dy_dx.subs({x:5}).n(10))#转浮点，总位数为10

输出：

4*t**3 + x/t
    3*x**2 + cos(x)
    cos(5) + 75
    75.28366219

至于高阶的微分，嗯，差点漏了

import sympy as sy
    x=sy.symbols('x')
    print(sy.diff(sy.sin(x),x,2))#二阶微分

输出：

-sin(x)

**二、积分**

积分就比较复杂了。有些函数是不能把积分写出来的。

import sympy as sy
    x=sy.symbols('x')#约定变量
    y=sy.sin(x)**2+4*x+3
    int_y_x=sy.integrate(y,x)#不定积分
    z=sy.sqrt(sy.exp(sy.sin(x))+2)
    int_z_x=sy.integrate(z,x)#这个不定积分写不出
    fix_int_y_x=sy.integrate(y,(x,0,1))#定积分
    #打印
    print(int_z_x)
    print(int_y_x)
    print(fix_int_y_x)
    print(fix_int_y_x.n())#默认转浮点的话是15位

输出：

Integral(sqrt(exp(sin(x)) + 2), x)
    2*x**2 + 7*x/2 - sin(x)*cos(x)/2
    -sin(1)*cos(1)/2 + 11/2
    5.27267564329358

不过有个好消息是，即使不定积分写不出来，但定积分的数值还是能写出来的，至于具体原因学过展开式的大家都懂。

import sympy as sy
    x=sy.symbols('x')#约定变量
    y=sy.exp(sy.sin(x)+sy.cos(x))
    int_y_x=sy.integrate(y,x)
    fix_int_y_x=sy.integrate(y,(x,0,1))
    print(int_y_x)
    print(fix_int_y_x)
    print(fix_int_y_x.n())

输出：

Integral(exp(sin(x))*exp(cos(x)), x)
    Integral(exp(sin(x))*exp(cos(x)), (x, 0, 1))
    3.70023523055316

当然了，有些人可能会想到用scipy求积分。

from scipy.integrate import quad,dblquad
    import numpy as np
    
    def f(x):
        return np.exp(np.sin(x)+np.cos(x))
    result,eps=quad(f,0,1)#定积分
    print(result)
    print(eps)#误差
    def f1(y,x):
        return 1
    def g(x):#y=-np.sqrt(1-x**2)是下边界函数
        return -np.sqrt(1-x**2)
    def h(x):#y=np.sqrt(1-x**2)是上边界函数
        return np.sqrt(1-x**2)
    res,ep=dblquad(f1,-1,1,g,h)#二重积分,-1<x<1
    print(res)#积分区域是个单位圆，积分函数是1的话，相当于求单位圆柱的体积
    print(np.pi)#有误差是正常的

输出：

3.70023523055316
    4.1080863494892926e-14
    3.141592653589797
    3.141592653589793

除了连续函数求积分，还能用辛普森法则和梯形法则来求，不过要提供列表或数组。

import numpy as np
    from scipy.integrate import quad,trapz,simps
    
    def f(x):
        return np.exp(np.sin(x)+np.cos(x))
    result,eps=quad(f,0,1)#定积分
    x=np.linspace(0,1,num=200)#200个点
    res=simps(f(x),x)#辛普森法则
    r=trapz(f(x),x)#梯形法则
    print(result)
    print(res)#点越多越精确
    print(r)

输出：

3.70023523055316
    3.700235203233112
    3.7002269868042608

在插值的基础上求积分也是个不错的应用。

更多花里胡哨的积分就要看官网了：

[https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html\#scipy.integrate.dblquad][https_docs.scipy.org_doc_scipy_reference_generated_scipy.integrate.dblquad.html_scipy.integrate.dblquad]

还有就是积分图像的描绘也是个不错的应用，除了plt还要引入Polygon函数。

from matplotlib import pyplot as plt
    import numpy as np
    from matplotlib.patches import Polygon
    
    def f(x):
        return np.exp(np.sin(x)+np.cos(x))
    x=np.linspace(-1,2,num=100)
    fig,ax=plt.subplots(figsize=(8,5))#绘制8*5的图形
    plt.plot(x,f(x))
    a,b=0,1#积分域0到1
    ix=np.linspace(a,b)
    iy=f(ix)
    verts=[(a,0)]+list(zip(ix,iy))+[(b,0)]#表示积分的区域
    #使用Polygon函数填充，表面颜色是0.7，边沿颜色是0.5
    poly=Polygon(verts,facecolor='0.7',edgecolor='0.5')
    ax.add_patch(poly)#将填充面积添加到图表里
    plt.text(0.5*(a+b),1,r'$\int_'+str(a)+'^'+str(b)+' f(x)\mathrm{d}x$',horizontalalignment='center',fontsize=20)
    #添加LaTex的标量名称需要用2个$符号包含在内才行,这标明是积分，水平居中对齐，坐标是（0.5，1），字体20
    plt.xlabel('$x$')
    plt.ylabel('$f(x)=exp(sin(x)*cos(x))$')
    ax.set_xticks((a,b))#x轴刻度位置
    ax.set_xticklabels(('$'+str(a)+'$','$'+str(b)+'$'))#内容
    ax.set_yticks([f(a),f(b)])
    ax.set_yticklabels(('$'+str(round(f(a),6))+'$','$'+str(round(f(b),6))+'$'))#刻度值四舍五入到小数点后6位
    plt.show()

输出：

![SouthEast][]

终于搞定，滚去睡觉。

[https_docs.scipy.org_doc_scipy_reference_generated_scipy.integrate.dblquad.html_scipy.integrate.dblquad]: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html#scipy.integrate.dblquad
[SouthEast]: /images/20220601/aa9abe16235547be9c4b254beddde7ba.png