Python在高等数学中的应用

写这篇博客的原因：
	在备考数学建模时，往往会遇见数学建模与高等数学相关的知识点结合的问题，这类问题用笔算是很简单的，但使用代码来实现求解就比较的麻烦啦。
	于是我就想整理一下在数学建模中以Python实现的代码中会用到那些工具库与相关语法。
下面我将分为5个方面介绍一下：

一、SymPy工具库的介绍

SymPy是Python版的开原计算机代数系统的实现，通俗的讲，SymPy是用于符号运算的工具库。现在这个工具库包括数十个模块，使用命令
>>>help("sympy")
可以看到这些模块的名称，具体的模块介绍

SymPy简介
符号运算的基础知识（较简略）：
使用该库进行符号计算，首先要建立符号变量以及符号表达式。
符号变量是构成符号表达式的基本元素，可通过库中的symbol()函数创建。
例如：

from sympy import *
x=symbols('x')    x是符号变量的名称，'x'是符号变量的值，用于显示
y,z=symbols('y z')    构建多个符号变量时，中间以空格分隔

下面写一段代码展示一下：

from sympy import *
x,y,z=symbols('x y z')
m0,m1,m2,m3=symbols('m0:4')  #创建多个符号变量
x=sin(1)
print("x=",x);print("x=",x.evalf())
print("x=",x.n(16))  # 显示小数点后16位数字
print("pi的两种显示格式:{},{}".format(pi,pi.evalf(3)))  # 这里不能使用n()函数
expr1=y*sin(y**2)  # 创建第一个符号表达式
expr2=y**2+sin(y)*cos(y)+sin(z)   # 创建第二个符号表达式
print("expr1=",expr1)
print("y=5时,expr1=",expr1.subs(y,5))  # 代入一个符号变量的值
print("y=2,z=3时，expr2=",expr2.subs({y:2,z:3}))
print("y=2,z=3时，expr2=",expr2.subs({y:2,z:3}).n())  # 以浮点数显示计算结果

运行结果为：

x= sin(1)
x= 0.841470984807897
x= 0.8414709848078965
pi的两种显示格式:pi,3.14
expr1= y*sin(y**2)
y=5时,expr1= 5*sin(25)
y=2,z=3时，expr2= sin(2)*cos(2) + sin(3) + 4
y=2,z=3时，expr2= 3.76271876040590

SymPy还有很多函数可以用于处理有理数。
下面介绍 together函数和apart函数相关操作：

from sympy import *
x1,x2,x3,x4=symbols('m1:5');x=symbols('x')
print(x1/x2+x3/x4)
print(together(x1/x2+x3/x4))
print((2*x**2+3*x+4)/(x+1))
print(simplify((2*x**2+3*x+4)/(x+1))) # 化简没有效果
#  请注意，在进行 SymPy 表达式化简时，并不能保证每个表达式都能被完全化简。某些复杂或特殊的表达式可能无法进一步简化。
print(apart((2*x**2+3*x+4)/(x+1)))   #  相当于分式的化简

运行结果：

m1/m2 + m3/m4
(m1*m4 + m2*m3)/(m2*m4)
(2*x**2 + 3*x + 4)/(x + 1)
(2*x**2 + 3*x + 4)/(x + 1)
2*x + 1 + 3/(x + 1)

关于该库的相关操作，需要自己不断的摸索，才可以掌握一些基础语法的。

二、SciPy工具库的简介

SciPy是对NumPy的功能扩展，其中涵盖了许多高级数学函数。
SciPy是在NumPy数组框架的基础上实现的，它对NumPy数组和基本的数组运算进行扩展，以满足相关的计算需求。
SciPy支持的功能包括文件处理、积分、数值分析、优化方法、统计学、信号与图像处理、聚类分析和空间分析等。

python之scipy库详解

三、用SumPy做符号函数画图

用SymPy做符号函数画图还是很方便的。
下面说明二维图形、三维图形和隐函数符号函数画图方法。

#  1、二维曲线画图
#  plot(表达式,变量取值范围,属性=属性值)
#  多重绘制的使用格式为
#  plot(表达式1,表达式2,变量取值范围,属性=属性值)
#  或者 plot((表达式1,变量取值范围1),(表达式2,变量取值范围2))

实例如下：

from sympy.plotting import plot
from sympy.abc import x,pi
from sympy.functions import sin,cos
plot((2*sin(x),(x,-6,6)),(cos(x+pi/4),(x,-5,5)),(sin(x**2),(x,-4,4)))

运行结果：

#  三维曲面画图
from pylab import rc
from sympy.plotting import plot3d
from sympy.abc import x,y
from sympy.functions import sin,sqrt
# rc('font',size=16);rc('text',usetex=True)
plot3d(sin(sqrt(x**2+y**2)),(x,-10,10),(y,-10,10),xlabel='$x$',ylabel="$y$")

运行结果：

#  3、隐函数画图    (x-1)^2+(y-2)^3-4=0
from pylab import  rc
from sympy import plot_implicit as pt,Eq
from sympy.abc import x,y
pt(Eq((x-1)**2+(y-2)**3,4),(x,-6,6),(y,-2,4),xlabel='$x$',
   ylabel='$y$')

运行结果：

对于最后一个，还可以使用匿名函数(lambda函数)设计

#  使用匿名函数设计
from sympy import plot_implicit as pt
from sympy.abc import x,y
ezplot=lambda  expr:pt(expr)
#这段代码中 `ezplot=lambda expr:pt(expr)` 是一个 lambda 函数，用于生成一个基于点集的散点图。
#该函数接受一个表达式 `expr` 作为输入，并使用 `pt()` 函数将其转换为点集。
# `pt()` 函数是 Matplotlib 库中的一个函数，用于绘制散点图。
#因此，这段代码的作用是将给定的表达式转换为散点图，并将其显示在图形界面上。
ezplot((x-1)**2+(y-2)**3-4)

运行的结果大致是一样的

四、高等数学问题的符号解

SymPy包括许多的功能。它主要处理三种类型的数据：整型数据、实数和有理数。

1、求极限

from sympy import *
x=symbols('x')
print(limit(sin(x)/x,x,0))
print(limit(pow(1+1/x,x),x,oo))  # 这里的  oo  表示正无穷

运行结果：

1
E

2、求导数

from sympy import *
x,y=symbols('x y')
z=sin(x)+x**2*exp(y)
# diff(func,x,n)
# 其中，func是要求导的函数，x是要对其求导的变量，n是可选的，表示求n阶导数，默认为1阶导数。
print("关于x的二阶偏导数为：",diff(z,x,2))
print("关于y的一阶偏导数为：",diff(z,y))
print("关于x的二阶偏导数,再关于y的二阶偏导数为：",diff(z,x,2,y,2))

运行结果：

关于x的二阶偏导数为： 2*exp(y) - sin(x)
关于y的一阶偏导数为： x**2*exp(y)
关于x的二阶偏导数,再关于y的二阶偏导数为： 2*exp(y)

将导数带入具体的值求某一点处的导数

from sympy import *
def func(x):
    return x**4
x = symbols("x")
print(diff(func(x),x).subs(x,2))   # 表示将x = 2，带入导函数4*x**3中

运行结果：

32

3、级数的求和

from sympy import *
k,n=symbols('k n')
print(summation(k**2,(k,1,n)))
print(factor(summation(k**2,(k,1,n))))
# 把计算结果因式分解
print(summation(1/k**2,(k,1,oo)))
# 这里是两个小“o”表示正无穷

运行结果：

n**3/3 + n**2/2 + n/6
n*(n + 1)*(2*n + 1)/6
pi**2/6

4、泰勒展开

from pylab import rc
from sympy import *
from sympy.plotting import *
# rc('font',size=16);rc('text',usetex=True)
x=symbols('x');y=sin(x)
# 定义了两个符号变量 x 和 y，其中 y 是一个三角函数的表达式。
for k in range(3,8,2):
# 使用 for 循环从 3 到 7（步长为 2）遍历每个整数 k。
    print(y.series(x,0,k)) # 等价于print(series(y,x,0,k))
# 在循环内部，使用 print 函数打印出 y 在给定的自变量 x 取值范围内的一系列近似表达式，直到取值达到 k。
plot(y,series(y,x,0,3).removeO(),series(y,x,0,5).removeO(),
     series(y,x,0,7).removeO(),(x,0,2),xlabel='$x$',ylabel="$y$")
#使用 plot 函数绘制 y 的曲线，其中包括使用 series 函数计算的从 k=0 到 k=3、k=5 和 k=7 的近似表达式，
# 以及在自变量 x 取值范围为 0 到 2 的区间内的数据点。
# 使用 xlabel 和 ylabel 函数设置坐标轴标签。

运行结果：

x + O(x**3)
x - x**3/6 + O(x**5)
x - x**3/6 + x**5/120 + O(x**7)

5、不定积分和定积分
python中Scipy模块求取积分的方法

from sympy import integrate,symbols,sin,pi,oo
x=symbols('x')
print(integrate(sin(2*x),(x,0,pi)))
print(integrate(sin(x)/x,(x,0,oo)))

运行结果：

0
pi/2

6、求解代数方程（组）的符号解

from sympy import *
x,y=symbols('x y')
print(solve(x**3-1,x))  #  解出来的解有虚数解
print(solve((x-2)**2*(x-1)**3,x))
print(roots((x-2)**2*(x-1)**3,x)) # roots可以得到根的重数信息

运行结果：

[1, -1/2 - sqrt(3)*I/2, -1/2 + sqrt(3)*I/2]
[1, 2]
{2: 2, 1: 3}

from sympy.abc import x,y
from sympy import solve
s=solve([x**2+y**2-1,x-y],[x,y])
print("方程组的解为: ",s)

运行结果：

方程组的解为:  [(-sqrt(2)/2, -sqrt(2)/2), (sqrt(2)/2, sqrt(2)/2)]

#  求函数f(x)=2x^3-5x^2+x的驻点，并求函数在[0,1]上的最大值
from sympy import *
x=symbols('x');y=2*x**3-5*x**2+x
x0=solve(diff(y,x),x) #求导后再求驻点
print("驻点的精确解为：",x0)
print("驻点的浮点数表示为：",[x0[i].n()for i in range(len(x0))])
# 列表中的符号数无法整体转换为浮点数
y0=[y.subs(x,0),y.subs(x,1),y.subs(x,x0[0]).n()]
# 带入区间端点和一个驻点的值
print("三个点的函数值分别为：",y0)

运行结果：

驻点的精确解为： [5/6 - sqrt(19)/6, sqrt(19)/6 + 5/6]
驻点的浮点数表示为： [0.106850176076554, 1.55981649059011]
三个点的函数值分别为： [0, -2, 0.0522051838383851]

7、求微分方程（方程组）的符号解
在SymPy库中有dsolve函数求常微分方程的符号解
在声明时，可以使用Function()函数
y=Function(‘y’) 或者 y=symbols(‘y’,cls=Function) 将符号变量声明为函数类型

# y"-5y'+6y=0;    y"-5y'+6y=xe^2x
from sympy import *
x=symbols('x');y=symbols('y',cls=Function)
eq1=diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)
eq2=diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)-x*exp(2*x)
print("齐次方程的解为：",dsolve(eq1,y(x)))
print("非齐次方程的解为：",dsolve(eq2,y(x)))

运行结果：

齐次方程的解为： Eq(y(x), (C1 + C2*exp(x))*exp(2*x))
非齐次方程的解为： Eq(y(x), (C1 + C2*exp(x) - x**2/2 - x)*exp(2*x))

#  求下列微分方程的解
#  y"-5y'+6y=0  y(0)=1,y'(0)=0
#  y"-5y'+6y=xe^2x  y(0)=1,y(2)=0
from sympy import *
x=symbols('x');y=symbols('y',cls=Function)  # 将符号变量y声明为函数类型
eq1=diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)
eq2=diff(y(x),x,2)-5*diff(y(x),x)+6*y(x)-x*exp(2*x)
print("初值问题的解为：{}".format(dsolve(eq1,y(x),ics={y(0):1,diff(y(x),x).subs(x,0):0})))
# 将初值条件y(0)=1和dy/dx|_(x=0)=0传递给了ics参数
y2=dsolve(eq2,y(x),ics={y(0):1,y(2):0})
print("边值问题的解为：{}".format(y2))

运行结果：

初值问题的解为：Eq(y(x), (3 - 2*exp(x))*exp(2*x))
边值问题的解为：Eq(y(x), (-x**2/2 - x + 3*exp(x)/(-1 + exp(2)) + (-4 + exp(2))/(-1 + exp(2)))*exp(2*x))

五、高等数学问题的数值解

大多数实际问题是无法求符号解的，只能求数值解，即近似解。
只介绍泰勒级数相关求解方法：

import numpy as np
import matplotlib.pyplot as plt
def fac(n): return (1 if n<1 else n*fac(n-1))
def item(n,x):return (-1)**n*x**(2*n+1)/fac(2*n+1)
def mysin(n,x):return (0 if n<0 else mysin(n-1,x)+item(n,x))
x=np.linspace(-2*np.pi,2*np.pi,101)
#   从-2π到2π之间取101个等距离的数值，然后用这些数值作为数组的元素值初始化x数组
plt.plot(x,np.sin(x),'*-')
#  '*-'：表示使用星号标记点并使用实线连接这些点，用于指定线条的样式。
str=['v-','H--','-.']
for n in [1,2,3]:plt.plot(x,mysin(2*n-1,x),str[n-1])
plt.legend(['sin','n=1','n=3','n=5'])
plt.savefig('figure3_16.png',dpi=500);plt.show()

运行结果：

总结

1、以上就是我在备考数学建模当中的获得的相关感悟，希望能帮到大家。
2、当然，上述的相关函数的功能肯定是不全的，需要大家备考时自己查阅资料。
3、博客中的疏漏在所难免，敬请大佬们不吝指正！！！