scipy 与 sympy 模块在数学建模中的应用

River Chandler

已于 2022-09-13 10:27:45 修改

阅读量1.2k

点赞数 5

分类专栏： # 数学模型与仿真文章标签： scipy

于 2022-09-11 00:25:17 首次发布

本文链接：https://blog.csdn.net/Chandler_river/article/details/126800889

版权

数学模型与仿真专栏收录该内容

10 篇文章 4 订阅

订阅专栏

也许会有人问你一个比较有趣的问题 sympy 是啥意思

import sympy

#Handbook of Symbolic Python

sympy可能很方便，但更有可能不能满足你的需求，如果只是想解决高数作业的话，没问题的。

1.定义变量

#define variables
x = sympy.Symbol("x")
y = sympy.Symbol("y")
z = sympy.Symbol("z")
#define variables end

2.定义函数

#define functions
f = x**2 + x**4 + x**2 + 2*x**4
g = x**3 + x * y**2 + x*y*z
h = x**2 + sympy.cos(y**2) + z**x
I = x*(x**2+x*y+x*y**2)
J = sympy.sin(x) + sympy.cos(y)
#define functions end

3.do it later 功能

#do it later
d = sympy.Derivative(h,x)
print(d)
print(d.doit())
#do it later end

4.多项式问题

因式分解；简化；展开；合并同类项；改写三角函数；泰勒展开；解一元方程

                            #polynomial arithmetic

#factoring
print("factor",end = " ")
print(sympy.factor(f))

#simplify
print("simplify",end= " ")
print(sympy.simplify(g))

#expand
print("expand",end = " ")
print(sympy.expand(I))

#collect
print("collect",end=" ")
print(sympy.collect(f,[x**2,x**4]))

#rewrite
print("rewrite",end = " ")
print(J.rewrite(sympy.tan(x)))

#series
print("series",end = " ")
print(J.series(x,0,10))

#solve
print("solve",end = " ")
print(x**2+x<10,x)

5.微积分

求极限；求导数和偏导数；积分；解常微分方程与常微分方程组

                            #calculus
#limit
print("limit",end = " ")
print(sympy.limit(sympy.sin(x)/x,x,0))

#calculate derivative and partial derivative
print(sympy.diff(g,z))
print(sympy.diff(f,x,3))

#integrate
print("sympy.integrate(f,x)",end = " ")
print(sympy.integrate(f,x))
print("sympy.integrate(f,(x,0,5))",end = " ")
print(sympy.integrate(f,(x,0,5)))

#dsolve
from sympy import Function,dsolve,Derivative,symbols,Eq
from sympy.abc import t
x = Function("x")
result = dsolve(Derivative(x(t),t,4)-22*Derivative(x(t),t,2)-24*x(t))
print("result",end="  ")
print(result)
print("ODEs")

x = Function("x")
y = Function("y")
eq=(Eq(Derivative(x(t),t,2),12*(x(t)+y(t))), \
    Eq(Derivative(y(t),t,2),12*x(t)+10*y(t)))
result = dsolve(eq)
print("ODEs-result",end=" ")
print(result)

6.矩阵运算

解线性方程组；解非线性方程组

                            #Matrix
from sympy.matrices import Matrix
x = symbols("x")
y = symbols("y")
result = sympy.linsolve([Eq(x-y,1),Eq(x+y,1)],(x,y))
print("linsolve",end=" ")
print(result)
result = sympy.nonlinsolve([Eq(x**2-y,1),Eq(x+y**2,3)],(x,y))
print("nonsolve",end=" ")
print(result)

7.所有代码的运行结果

Derivative(x**2 + z**x + cos(y**2), x)
2*x + z**x*log(z)
factor x**2*(3*x**2 + 2)
simplify x*(x**2 + y**2 + y*z)
expand x**3 + x**2*y**2 + x**2*y
collect 3*x**4 + 2*x**2
rewrite (1 - tan(y/2)**2)/(tan(y/2)**2 + 1) + 2*tan(x/2)/(tan(x/2)**2 + 1)
series cos(y) + x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)
solve x**2 + x < 10 x
limit 1
x*y
72*x
sympy.integrate(f,x) 3*x**5/5 + 2*x**3/3
sympy.integrate(f,(x,0,5)) 5875/3
result  Eq(x(t), C1*exp(-t*sqrt(11 + sqrt(145))) + C2*exp(t*sqrt(11 + sqrt(145))) + C3*sin(t*sqrt(-11 + sqrt(145))) + C4*cos(t*sqrt(-11 + sqrt(145))))
ODEs
ODEs-result [Eq(x(t), C1*sqrt(11 + sqrt(145))*(67 - 5*sqrt(145))*exp(t*sqrt(11 + sqrt(145)))/144 - C2*sqrt(-11 + sqrt(145))*(5*sqrt(145) + 67)*sin(t*sqrt(-11 + sqrt(145)))/144 - C3*sqrt(-11 + sqrt(145))*(5*sqrt(145) + 67)*cos(t*sqrt(-11 + sqrt(145)))/144 - C4*sqrt(11 + sqrt(145))*(67 - 5*sqrt(145))*exp(-t*sqrt(11 + sqrt(145)))/144), Eq(y(t), -C1*(11 - sqrt(145))*sqrt(11 + sqrt(145))*exp(t*sqrt(11 + sqrt(145)))/24 + C2*sqrt(-11 + sqrt(145))*(11 + sqrt(145))*sin(t*sqrt(-11 + sqrt(145)))/24 + C3*sqrt(-11 + sqrt(145))*(11 + sqrt(145))*cos(t*sqrt(-11 + sqrt(145)))/24 + C4*(11 - sqrt(145))*sqrt(11 + sqrt(145))*exp(-t*sqrt(11 + sqrt(145)))/24)]
linsolve {(1, 0)}
nonsolve {(3 - (-sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2 - sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) - 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2)**2, -sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2 - sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) - 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2), (3 - (-sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2 + sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) - 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2)**2, -sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2 + sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) - 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2), (3 - (-sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2 + sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2)**2, -sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2 + sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2), (3 - (sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2 + sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2)**2, sqrt(-2*(sqrt(12081)/144 + 113/16)**(1/3) - 22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2/sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4) + 8)/2 + sqrt(22/(3*(sqrt(12081)/144 + 113/16)**(1/3)) + 2*(sqrt(12081)/144 + 113/16)**(1/3) + 4)/2)}

8.subs()函数

from sympy import symbols,diff

def fun(x,y):
    return x**2 + 2*y**2

x,y = symbols("x y",real=True)
print("(partial f/partial x)",end="  ")
print(diff(fun(x,y),x))
print("numeric (partial f/partial x)",end="  ")
print(diff(fun(x,y),x).subs({x:3,y:1}))

(partial f/partial x)  2*x
numeric (partial f/partial x)  6

scipy 指的是Science Python

1.scipy 做积分与重积分

import scipy.integrate as integrate
import numpy as np

#一重积分的lambda表达法
result1 = integrate.quad(lambda x:x**2+np.exp(x),0,5)
print("result1",end=" ")
print(result1)


#一重积分的函数表达法
a,b,c = 10,5,5
ARG = (a,b,c)
def fun(x,a,b,c):
    return a*x*(x>5)+b*x*(x<4)+c*x*(4<=x<=5)

result2 = integrate.quad(fun,3,6,args=ARG)
print("result2",end=" ")
print(result2)


#多重积分
def fun2(x,y,z,t,a,b):
    return x*y*z*t*a+b

def bounds_x(*args):
    return [0,3]
def bounds_y(*args):
    return [0,3]
def bounds_z(*args):
    return [0,0.5]
def bounds_t(*args):
    return [0,3]

a,b = 5,5
ARG = (a,b)

result3 = integrate.nquad(fun2,[bounds_x,bounds_y,bounds_z,bounds_t],args=ARG)
print("result3",end=" ")
print(result3)

运行结果

result1 (189.07982576924329, 2.099207760611269e-12)
result2 (95.00000000000001, 1.0547118733938988e-13)
result3 (124.45312500000003, 2.1784368501755553e-12)

2.scipy 计算行列式与逆矩阵

from scipy import linalg
import numpy as np

arr = np.array([[1,2],
               [3,4]])
#det
print("det",end=" ")
result = linalg.det(arr)
print(result)
#det end

#inv
result = linalg.inv(arr)
print("inv",end=" ")
print(result)
#inv end

det -2.0
inv [[-2.   1. ]
 [ 1.5 -0.5]]

3.scipy 单变量插值

#interpolate
#单变量插值UnivariateSpline
from scipy import interpolate
#x,y是X-Y坐标数组
#w是每个数据点的权重值
#k为样条曲线的阶数
#s为平滑参数

x = np.linspace(0,10,20)
y = np.sin(x)

sx = np.linspace(0,12,100)
func1 = interpolate.UnivariateSpline(x,y,s=0)
func2 = interpolate.UnivariateSpline(x,y)
sy1 = func1(sx)
sy2 = func2(sx)
import matplotlib.pyplot as plt
#plt.plot(x,y,"o")
#plt.plot(sx,sy1)
#plt.plot(sx,sy2)
#plt.show()

4.scipy 二维插值

#二维插值Rbf()
def func(x,y):
    return (x+y)*np.exp(-5*(x**2+y**2))

x = np.random.uniform(low=-1,high=1,size=100)
y = np.random.uniform(low=-1,high=1,size=100)
z = func(x,y)
func=interpolate.Rbf(x,y,z,function='multiquadric')
xnew,ynew=np.mgrid[-1:1:100j,-1:1:100j]
znew=func(xnew,ynew)#输入输出都是二维

import mpl_toolkits.mplot3d
import matplotlib.pyplot as plt
ax=plt.subplot(111,projection='3d')
ax.plot_surface(xnew,ynew,znew)
ax.scatter(x,y,z,c='r',marker='^')
plt.savefig("0.jpg")

5.scipy solve_ivp 解一类特殊的一阶微分方程和一阶微分方程组
- $\frac{dy(x)}{dx}=...$
- $\begin{Bmatrix} \frac{dy_1(x)}{dx}=...\\ \frac{df_2(x)}{dx}=...\\ ... \end{Bmatrix}$

def fun(t,y):
    y1 = y[0]
    y2 = y[1]
    return [y2+t,y1*y2-t]

# 初始条件
y0 = [0,1]
y = solve_ivp(fun,(0,50),y0,method='Radau',t_eval=np.arange(1,50,1))
t = y.t
data = y.y
plt.plot(t,data[0],label="y1",c="blue")
plt.plot(t,data[1],label="y2",c="red")

plt.xlabel("t")
plt.legend()
plt.savefig("0.jpg")

def fun(t,y):
    return [t+y[0]]
    
# 初始条件
y0 = [0]
y = solve_ivp(fun,(0,50),y0,method='Radau',t_eval=np.arange(1,50,1))
t = y.t
data = y.y
plt.plot(t, data[0])
plt.xlabel("t")
plt.legend()
plt.savefig("0.jpg")

6. scipy.misc.derivative 求单变量函数的导数

from scipy.misc import derivative
def f(x):
    return x**2+x**3+x
def df(x):
    return derivative(f,x)
print(f(1))
print(df(1))

运行结果

3
7.0

7.scipy curve_fit 拟合
- 拟合和插值的区别
  - 插值
    - 只管搞出一个平面。不需要管平面什么样
  - 拟合
    - 搞出一个基本有框架的平面

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

#确定你想要的函数
def func(x,a,b,c):
    return a * x**2 * np.exp(b*x) + c * x**0.5
#end

#生产数据
x_data = np.linspace(0,10,500)
y_data = func(x_data,2.5,1.3,0.5)
np.random.seed(0)
err_stdev = 0.2
y_noise = err_stdev * np.random.normal(size=x_data.size)
y_data += y_noise
#end

#绘制原图
plt.title("curve_fit")
plt.plot(x_data,y_data,"r-.",label="raw data")
#end


#无范围拟合
popt,pcov = curve_fit(func,x_data,y_data)
plt.plot(x_data,func(x_data,*popt),"b--",label="fit first")
print("popt 1",end=" ")
print(popt)
print("pcov 1")
print(pcov)
#end

#限定参数范围拟合
popt,pcov = curve_fit(func,x_data,y_data,bounds=([0,0.2,0.4],[4,0.4,0.8]))
plt.plot(x_data,func(x_data,*popt),"y-",label="fit second")
print("popt 2",end=" ")
print(popt)
print("pcov 2")
print(pcov)
#end

plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

#end

协方差pcov表示的是两个变量总体误差的期望

popt 1 [2.49999999 1.3        0.49486601]
pcov 1
[[ 1.44050764e-15 -5.88808815e-17 -1.06696225e-10]
 [-5.88808815e-17  2.40940547e-18  4.24978186e-12]
 [-1.06696225e-10  4.24978186e-12  2.84981426e-05]]
popt 2 [4.  0.4 0.8]
pcov 2
[[ 1.29983692e+08 -3.29981981e+06 -7.68847399e+09]
 [-3.29981981e+06  8.42235301e+04  1.87368753e+08]
 [-7.68847399e+09  1.87368753e+08  7.57765000e+11]]

多变量拟合

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

#确定你想要的函数
def func(X,a,b,c):
    x,y = X
    return a * x**2 * np.exp(b*x) + c * x**0.5 * y
#end

#生产数据
x_data = np.linspace(0,10,100)
y_data = np.linspace(0,10,100)
z_data = func((x_data,y_data),2.5,1.3,0.5)
np.random.seed(0)
err_stdev = 0.2
z_noise = err_stdev * np.random.normal(size=x_data.size)
z_data += z_noise
#end


popt,pcov = curve_fit(func,(x_data,y_data),z_data)
print("popt  ")
print(popt)
print("pcov  ")
print(pcov)

输出

popt  
[2.50000018 1.29999999 0.49521289]
pcov  
[[ 8.95177268e-15 -3.60828062e-16 -1.48495956e-10]
 [-3.60828062e-16  1.45592255e-17  5.85915343e-12]
 [-1.48495956e-10  5.85915343e-12  5.05947283e-06]]

8.optimize.minimize 有约束，无约束寻找最小值

from scipy.optimize import minimize
import numpy as np

#无约束寻找最小值
def f(x):
    
    return x[0]**2+(x[1]-3)*(x[2]+2)+x[0]*2+x[1]*0.5+x[2]**2
x0 = [4,2,3]

result = minimize(f,x0)
print(result)

print("###############")
#有约束寻找最小值
cons = ({'type': 'eq',
         'fun': lambda x: np.array([x[0] ** 3 - x[1]*x[2]])},
        {'type': 'ineq',
         'fun': lambda x: np.array([x[1] - 1])})
result = minimize(f,x0,constraints=cons,method="SLSQP")
print(result)

运行结果

      fun: -3797921.0484094834
 hess_inv: array([[ 3.03316132,  5.23929226, -1.34256663],
       [ 5.23929226,  9.70645543, -1.51875968],
       [-1.34256663, -1.51875968,  0.78806632]])
      jac: array([-9175.  ,  2385.5 , -8028.25])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 472
      nit: 4
     njev: 115
   status: 2
  success: False
        x: array([ -4588.41338386, -12791.39118934,   2382.99401134])
###############
     fun: -4.009968593644863
     jac: array([ 1.15420681,  2.42436856, -2.15126294])
 message: 'Optimization terminated successfully'
    nfev: 92
     nit: 21
    njev: 21
  status: 0
 success: True
       x: array([-0.4228966 ,  1.        , -0.07563147])

9. root 寻找一元函数的根

from scipy.optimize import root
import numpy as np
print(root(lambda x:x*(x-2)-1,4))

    fjac: array([[-1.]])
     fun: array([5.99520433e-15])
 message: 'The solution converged.'
    nfev: 9
     qtf: array([-3.68390629e-09])
       r: array([-2.82843182])
  status: 1
 success: True
       x: array([2.41421356])