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

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

    1. import sympy
    2.  
    3. #Handbook of Symbolic Python

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

    • 1.定义变量

    1. #define variables
    2. x = sympy.Symbol("x")
    3. y = sympy.Symbol("y")
    4. z = sympy.Symbol("z")
    5. #define variables end

    • 2.定义函数

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

    • 3.do it later 功能

    1. #do it later
    2. d = sympy.Derivative(h,x)
    3. print(d)
    4. print(d.doit())
    5. #do it later end

    • 4.多项式问题

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

    1.                             #polynomial arithmetic
    2.  
    3. #factoring
    4. print("factor",end = " ")
    5. print(sympy.factor(f))
    6.  
    7. #simplify
    8. print("simplify",end= " ")
    9. print(sympy.simplify(g))
    10.  
    11. #expand
    12. print("expand",end = " ")
    13. print(sympy.expand(I))
    14.  
    15. #collect
    16. print("collect",end=" ")
    17. print(sympy.collect(f,[x**2,x**4]))
    18.  
    19. #rewrite
    20. print("rewrite",end = " ")
    21. print(J.rewrite(sympy.tan(x)))
    22.  
    23. #series
    24. print("series",end = " ")
    25. print(J.series(x,0,10))
    26.  
    27. #solve
    28. print("solve",end = " ")
    29. print(x**2+x<10,x)

    • 5.微积分

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

    1.                             #calculus
    2. #limit
    3. print("limit",end = " ")
    4. print(sympy.limit(sympy.sin(x)/x,x,0))
    5.  
    6. #calculate derivative and partial derivative
    7. print(sympy.diff(g,z))
    8. print(sympy.diff(f,x,3))
    9.  
    10. #integrate
    11. print("sympy.integrate(f,x)",end = " ")
    12. print(sympy.integrate(f,x))
    13. print("sympy.integrate(f,(x,0,5))",end = " ")
    14. print(sympy.integrate(f,(x,0,5)))
    15.  
    16. #dsolve
    17. from sympy import Function,dsolve,Derivative,symbols,Eq
    18. from sympy.abc import t
    19. x = Function("x")
    20. result = dsolve(Derivative(x(t),t,4)-22*Derivative(x(t),t,2)-24*x(t))
    21. print("result",end="  ")
    22. print(result)
    23. print("ODEs")
    24.  
    25. x = Function("x")
    26. y = Function("y")
    27. eq=(Eq(Derivative(x(t),t,2),12*(x(t)+y(t))), \
    28.     Eq(Derivative(y(t),t,2),12*x(t)+10*y(t)))
    29. result = dsolve(eq)
    30. print("ODEs-result",end=" ")
    31. print(result)

    • 6.矩阵运算

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

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

    • 7.所有代码的运行结果

    1. Derivative(x**2 + z**x + cos(y**2), x)
    2. 2*x + z**x*log(z)
    3. factor x**2*(3*x**2 + 2)
    4. simplify x*(x**2 + y**2 + y*z)
    5. expand x**3 + x**2*y**2 + x**2*y
    6. collect 3*x**4 + 2*x**2
    7. rewrite (1 - tan(y/2)**2)/(tan(y/2)**2 + 1) + 2*tan(x/2)/(tan(x/2)**2 + 1)
    8. series cos(y) + x - x**3/6 + x**5/120 - x**7/5040 + x**9/362880 + O(x**10)
    9. solve x**2 + x < 10 x
    10. limit 1
    11. x*y
    12. 72*x
    13. sympy.integrate(f,x) 3*x**5/5 + 2*x**3/3
    14. sympy.integrate(f,(x,0,5)) 5875/3
    15. 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))))
    16. ODEs
    17. 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)]
    18. linsolve {(1, 0)}
    19. 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()函数 

    1. from sympy import symbols,diff
    2.  
    3. def fun(x,y):
    4.     return x**2 + 2*y**2
    5.  
    6. x,y = symbols("x y",real=True)
    7. print("(partial f/partial x)",end="  ")
    8. print(diff(fun(x,y),x))
    9. print("numeric (partial f/partial x)",end="  ")
    10. print(diff(fun(x,y),x).subs({x:3,y:1}))
    1. (partial f/partial x)  2*x
    2. numeric (partial f/partial x)  6

    scipy  指的是Science Python

    • 1.scipy 做积分与重积分

    1. import scipy.integrate as integrate
    2. import numpy as np
    3.  
    4. #一重积分的lambda表达法
    5. result1 = integrate.quad(lambda x:x**2+np.exp(x),0,5)
    6. print("result1",end=" ")
    7. print(result1)
    8.  
    9.  
    10. #一重积分的函数表达法
    11. a,b,c = 10,5,5
    12. ARG = (a,b,c)
    13. def fun(x,a,b,c):
    14.     return a*x*(x>5)+b*x*(x<4)+c*x*(4<=x<=5)
    15.  
    16. result2 = integrate.quad(fun,3,6,args=ARG)
    17. print("result2",end=" ")
    18. print(result2)
    19.  
    20.  
    21. #多重积分
    22. def fun2(x,y,z,t,a,b):
    23.     return x*y*z*t*a+b
    24.  
    25. def bounds_x(*args):
    26.     return [0,3]
    27. def bounds_y(*args):
    28.     return [0,3]
    29. def bounds_z(*args):
    30.     return [0,0.5]
    31. def bounds_t(*args):
    32.     return [0,3]
    33.  
    34. a,b = 5,5
    35. ARG = (a,b)
    36.  
    37. result3 = integrate.nquad(fun2,[bounds_x,bounds_y,bounds_z,bounds_t],args=ARG)
    38. print("result3",end=" ")
    39. print(result3)

    运行结果

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

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

    1. from scipy import linalg
    2. import numpy as np
    3.  
    4. arr = np.array([[1,2],
    5.                [3,4]])
    6. #det
    7. print("det",end=" ")
    8. result = linalg.det(arr)
    9. print(result)
    10. #det end
    11.  
    12. #inv
    13. result = linalg.inv(arr)
    14. print("inv",end=" ")
    15. print(result)
    16. #inv end
    1. det -2.0
    2. inv [[-2.   1. ]
    3.  [ 1.5 -0.5]]

    • 3.scipy 单变量插值

    1. #interpolate
    2. #单变量插值UnivariateSpline
    3. from scipy import interpolate
    4. #x,y是X-Y坐标数组
    5. #w是每个数据点的权重值
    6. #k为样条曲线的阶数
    7. #s为平滑参数
    8.  
    9. x = np.linspace(0,10,20)
    10. y = np.sin(x)
    11.  
    12. sx = np.linspace(0,12,100)
    13. func1 = interpolate.UnivariateSpline(x,y,s=0)
    14. func2 = interpolate.UnivariateSpline(x,y)
    15. sy1 = func1(sx)
    16. sy2 = func2(sx)
    17. import matplotlib.pyplot as plt
    18. #plt.plot(x,y,"o")
    19. #plt.plot(sx,sy1)
    20. #plt.plot(sx,sy2)
    21. #plt.show()

    • 4.scipy  二维插值

    1. #二维插值Rbf()
    2. def func(x,y):
    3.     return (x+y)*np.exp(-5*(x**2+y**2))
    4.  
    5. x = np.random.uniform(low=-1,high=1,size=100)
    6. y = np.random.uniform(low=-1,high=1,size=100)
    7. z = func(x,y)
    8. func=interpolate.Rbf(x,y,z,function='multiquadric')
    9. xnew,ynew=np.mgrid[-1:1:100j,-1:1:100j]
    10. znew=func(xnew,ynew)#输入输出都是二维
    11.  
    12. import mpl_toolkits.mplot3d
    13. import matplotlib.pyplot as plt
    14. ax=plt.subplot(111,projection='3d')
    15. ax.plot_surface(xnew,ynew,znew)
    16. ax.scatter(x,y,z,c='r',marker='^')
    17. 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}
    1. def fun(t,y):
    2.     y1 = y[0]
    3.     y2 = y[1]
    4.     return [y2+t,y1*y2-t]
    5.  
    6. # 初始条件
    7. y0 = [0,1]
    8. y = solve_ivp(fun,(0,50),y0,method='Radau',t_eval=np.arange(1,50,1))
    9. t = y.t
    10. data = y.y
    11. plt.plot(t,data[0],label="y1",c="blue")
    12. plt.plot(t,data[1],label="y2",c="red")
    13.  
    14. plt.xlabel("t")
    15. plt.legend()
    16. plt.savefig("0.jpg")

    1. def fun(t,y):
    2.     return [t+y[0]]
    3.     
    4. # 初始条件
    5. y0 = [0]
    6. y = solve_ivp(fun,(0,50),y0,method='Radau',t_eval=np.arange(1,50,1))
    7. t = y.t
    8. data = y.y
    9. plt.plot(t, data[0])
    10. plt.xlabel("t")
    11. plt.legend()
    12. plt.savefig("0.jpg")

     

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

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

    运行结果

    1. 3
    2. 7.0

    • 7.scipy curve_fit 拟合

      • 拟合和插值的区别
        • 插值
          • 只管搞出一个平面。不需要管平面什么样
        • 拟合
          • 搞出一个基本有框架的平面
    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy.optimize import curve_fit
    4.  
    5. #确定你想要的函数
    6. def func(x,a,b,c):
    7.     return a * x**2 * np.exp(b*x) + c * x**0.5
    8. #end
    9.  
    10. #生产数据
    11. x_data = np.linspace(0,10,500)
    12. y_data = func(x_data,2.5,1.3,0.5)
    13. np.random.seed(0)
    14. err_stdev = 0.2
    15. y_noise = err_stdev * np.random.normal(size=x_data.size)
    16. y_data += y_noise
    17. #end
    18.  
    19. #绘制原图
    20. plt.title("curve_fit")
    21. plt.plot(x_data,y_data,"r-.",label="raw data")
    22. #end
    23.  
    24.  
    25. #无范围拟合
    26. popt,pcov = curve_fit(func,x_data,y_data)
    27. plt.plot(x_data,func(x_data,*popt),"b--",label="fit first")
    28. print("popt 1",end=" ")
    29. print(popt)
    30. print("pcov 1")
    31. print(pcov)
    32. #end
    33.  
    34. #限定参数范围拟合
    35. popt,pcov = curve_fit(func,x_data,y_data,bounds=([0,0.2,0.4],[4,0.4,0.8]))
    36. plt.plot(x_data,func(x_data,*popt),"y-",label="fit second")
    37. print("popt 2",end=" ")
    38. print(popt)
    39. print("pcov 2")
    40. print(pcov)
    41. #end
    42.  
    43. plt.xlabel("x")
    44. plt.ylabel("y")
    45. plt.legend()
    46. plt.show()
    47.  
    48. #end

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

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

     多变量拟合

    1. import matplotlib.pyplot as plt
    2. import numpy as np
    3. from scipy.optimize import curve_fit
    4.  
    5. #确定你想要的函数
    6. def func(X,a,b,c):
    7.     x,y = X
    8.     return a * x**2 * np.exp(b*x) + c * x**0.5 * y
    9. #end
    10.  
    11. #生产数据
    12. x_data = np.linspace(0,10,100)
    13. y_data = np.linspace(0,10,100)
    14. z_data = func((x_data,y_data),2.5,1.3,0.5)
    15. np.random.seed(0)
    16. err_stdev = 0.2
    17. z_noise = err_stdev * np.random.normal(size=x_data.size)
    18. z_data += z_noise
    19. #end
    20.  
    21.  
    22. popt,pcov = curve_fit(func,(x_data,y_data),z_data)
    23. print("popt  ")
    24. print(popt)
    25. print("pcov  ")
    26. print(pcov)

    输出

    1. popt  
    2. [2.50000018 1.29999999 0.49521289]
    3. pcov  
    4. [[ 8.95177268e-15 -3.60828062e-16 -1.48495956e-10]
    5.  [-3.60828062e-16  1.45592255e-17  5.85915343e-12]
    6.  [-1.48495956e-10  5.85915343e-12  5.05947283e-06]]

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

    1. from scipy.optimize import minimize
    2. import numpy as np
    3.  
    4. #无约束寻找最小值
    5. def f(x):
    6.     
    7.     return x[0]**2+(x[1]-3)*(x[2]+2)+x[0]*2+x[1]*0.5+x[2]**2
    8. x0 = [4,2,3]
    9.  
    10. result = minimize(f,x0)
    11. print(result)
    12.  
    13. print("###############")
    14. #有约束寻找最小值
    15. cons = ({'type': 'eq',
    16.          'fun': lambda x: np.array([x[0] ** 3 - x[1]*x[2]])},
    17.         {'type': 'ineq',
    18.          'fun': lambda x: np.array([x[1] - 1])})
    19. result = minimize(f,x0,constraints=cons,method="SLSQP")
    20. print(result)

    • 运行结果

    1.       fun: -3797921.0484094834
    2.  hess_inv: array([[ 3.03316132,  5.23929226, -1.34256663],
    3.        [ 5.23929226,  9.70645543, -1.51875968],
    4.        [-1.34256663, -1.51875968,  0.78806632]])
    5.       jac: array([-9175.  ,  2385.5 , -8028.25])
    6.   message: 'Desired error not necessarily achieved due to precision loss.'
    7.      nfev: 472
    8.       nit: 4
    9.      njev: 115
    10.    status: 2
    11.   success: False
    12.         x: array([ -4588.41338386, -12791.39118934,   2382.99401134])
    13. ###############
    14.      fun: -4.009968593644863
    15.      jac: array([ 1.15420681,  2.42436856, -2.15126294])
    16.  message: 'Optimization terminated successfully'
    17.     nfev: 92
    18.      nit: 21
    19.     njev: 21
    20.   status: 0
    21.  success: True
    22.        x: array([-0.4228966 ,  1.        , -0.07563147])

    • 9. root 寻找一元函数的根

    1. from scipy.optimize import root
    2. import numpy as np
    3. print(root(lambda x:x*(x-2)-1,4))
    4.  
    5.     fjac: array([[-1.]])
    6.      fun: array([5.99520433e-15])
    7.  message: 'The solution converged.'
    8.     nfev: 9
    9.      qtf: array([-3.68390629e-09])
    10.        r: array([-2.82843182])
    11.   status: 1
    12.  success: True
    13.        x: array([2.41421356])

  • 相关阅读:
    【Tomcat专题】如何正确使用线程池
    神经网络和pid有什么区别,基于神经网络的pid控制
    Golang1.21更新内容全面介绍~
    深圳市商务局2022年度中央资金（跨境电子商务企业市场开拓扶持事项）申报指南
    17.4、JavaWeb-HTML、W3C、基本标签、html的各种标签
    【小航的算法日记】图论
    golang的垃圾回收算法之五GMP模型
    2023系统架构师---论软件系统架构评估（范文）
    设计模式-代理模式-笔记
    BloomFilter 布隆过滤器
  • 原文地址：https://blog.csdn.net/Chandler_river/article/details/126800889