In [1]:
import numpy as np
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
def fun(t,y):
"""System of ODEs."""
return -2*y +1
t_span = (0, 5) # time span
y0 = [0] # initial condition
sol = solve_ivp(fun, t_span, y0, method='RK45', t_eval=np.linspace(0, 5, 100))
#显示结果
plt.figure(figsize=(8, 4))
plt.plot(sol.t, sol.y[0],'r--', label='y(t)')
plt.title('Solution of ODE using solve_ivp')
plt.xlabel('t')
plt.ylabel('y')
plt.show()
例2 求解以下常微分方程 $y^{'''} - 3*y^{''} - y^{'}*y = 0$ 其边界值条件为$y^{'}(0)=1,y^{''}(0)=-1$
解：
- 设$y_{1}=y,y_{2} = y^{'},y_{3}=y^{'''}$那么常微分方程变成以下方程组: $$ \left\{ \begin{aligned} y_{1}^{'} = y_{2} \\ y_{2}^{'} = y_{3} \\ y_{3}^{'} = 3*y_{3}+ y_{2}*y_{1}\\ \end{aligned} \right. $$
In [2]:
import numpy as np
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
def fun2(x,Y):
y1,y2,y3 = Y # u1,w1分别表示y1,y2
dy1 = y2
dy2 = y3
dy3 = 3*y3 + y2*y1
return [dy1, dy2,dy3]
Y0 = [0,1,-1]
# 初始条件
x_span = (0, 1) # 时间区间
x_eval = np.linspace(0, 1, 100) # 计算点
sol = solve_ivp(fun2, x_span, Y0, method='RK45', t_eval=x_eval)
# 显示结果
plt.plot(sol.t, sol.y[0],'r-',label='y(t)')
plt.plot(sol.t, sol.y[1], 'g--',label="y'(t)")
plt.plot(sol.t, sol.y[2], 'b-.',label="y''(t)")
plt.legend()
plt.xlabel('t')
plt.show()
In [4]:
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
# 一维插值
x = np.linspace(0, 10, 10)
y = np.sin(x)
f_linear = interp1d(x, y)
f_cubic = interp1d(x, y, kind='cubic')
xnew = np.linspace(0, 10, 100)
plt.plot(x, y, 'o', label='data')
plt.plot(xnew, f_linear(xnew), 'r-', label='linear')
plt.plot(xnew, f_cubic(xnew), 'b--', label='cubic')
plt.legend()
plt.show()
- 曲线拟合
In [5]:
from scipy.optimize import curve_fit
from matplotlib import pyplot as plt
import numpy as np
def func(x, a, b, c):
return a * np.exp(b * x) + c
xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 2) + 2 * np.random.normal(size=len(xdata))
popt, pcov = curve_fit(func, xdata, y)
# print('Fitted parameters:', popt)
plt.plot(xdata, y, 'rp', label='benchmark data')
plt.plot(xdata, func(xdata,popt[0],popt[1],popt[2]), 'b--', label='curve_fit')
plt.legend()
plt.show()
3 线性规划¶
scipy.optimize.linprog是求解线性规划问题的工具- 例: 线性规划 $$ minz = -1*x_{1} - 2*x_{2} \\ s.t. \left\{\begin{aligned} x_{1} +2 *x_{2} \leq 4 \\ 4*x_{1} + 2*x_{2} \leq 12 \\ x_{1},x_{2} \geq 0\\ \end{aligned}\right. $$
In [8]:
from scipy.optimize import linprog
# 最小化 -x - 2y (即最大化 x + 2y)
c = [-1, -2]
# 不等式约束: x + 2y <= 4, 4x + 2y <= 12
A = [[1, 2],
[4, 2]]
b = [4, 12]
# 变量 x, y >=0
res = linprog(c, A_ub=A, b_ub=b, bounds=[(0, None), (0, None)])
print('Optimal value:', -res.fun) # 因为linprog默认是minimize
print('x:', res.x[0])
print('y:', res.x[1])
Optimal value: 4.0 x: 0.0 y: 2.0
4 非线性规划¶
scipy.optimize.minimize支持多种方法求解非线性优化。 $$ minz = 0.2*(x_1^2+x_2^2+x_3^2)+5*(x_1+x_2+x_3)+4*(2*x_1+x_2-140) \\ \left \{ \begin{array}{ll} x_1\geq 40 \\x_1+x_2 \geq 100 \\ x_1 +x_2+x_3 \geq 180 \\ 0 \leq x_1,x_2,x_3 \leq 100 \end{array} \right. $$
In [12]:
from scipy.optimize import minimize
def fun(x):
return 0.2*(x[0]**2+x[1]**2+x[2]**2)+5*(x[0]+x[1]+x[2])+4*(2*x[0]+x[1]-140)
x0 = [40, 60,80] # 初始猜测
cons= ({'type':'ineq', 'fun': lambda x: -100 + x[0]+ x[1]},
{'type':'ineq', 'fun': lambda x:- 180+ x[0]+ x[1]+x[2]})
# 约束条件
bounds = ((40, 100), (0, 100), (0, 100))
res = minimize(fun, x0, constraints=cons,bounds=bounds)
print('Optimal value:{:2.0f}'.format(res.fun) )
print('Optimal x: {:3.0f},{:3.0f},{:3.0f}'.format(res.x[0],res.x[1],res.x[2]))
Optimal value:3180 Optimal x: 50, 60, 70