python之微积分

一、微分

微分的话其实比较简单，手工算其实可以的，但是遇到特别复杂的函数，还是计算机算比较省时间，用到了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

还有就是积分图像的描绘也是个不错的应用，除了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()

输出：

终于搞定，滚去睡觉。