- 求下列微分方程的数值解,画出解的图形。
- spyder、matlab代码
- (1)spyder
-
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
def dy_dt(y,t): # 定义函数 f(y,t)
return np.sin(t**2)
y0 = [1] # y0 = 1 也可以
t=np.arange(-10,10,0.01) # (start,stop,step)
y=odeint(dy_dt,y0,t) # 求解微分方程初值问题
# 绘图
plt.plot(t,y)
plt.title('scipy.integrate.odeint')
plt.show()
print(y)
from scipy.integrate import odeint # 导入 scipy.integrate 模块 import numpy as np import matplotlib.pyplot as plt def dy_dt(y,t): # 定义函数 f(y,t) return np.sin(t**2) y0 = [1] # y0 = 1 也可以 t=np.arange(-10,10,0.01) # (start,stop,step) y=odeint(dy_dt,y0,t) # 求解微分方程初值问题 # 绘图 plt.plot(t,y) plt.title('scipy.integrate.odeint') plt.show() print(y)
-
(2)spyder
-
from scipy.integrate import odeint # 导入 scipy.integrate 模块
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 导数函数, 求 W=[x,y,z] 点的导数 dW/dt
def lorenz(W,t): # by youcans
p,r,b=10,28,3
x, y, z = W # W=[x,y,z]
dx_dt = p*(y-x) # dx/dt = p*(y-x), p: sigma
dy_dt = x*(r-z)-y # dy/dt = x*(r-z)-y, r:rho
dz_dt = x*y - b*z # dz/dt = x*y - b*z, b;beta
return np.array([dx_dt, dy_dt, dz_dt]);
t = np.arange(0, 30, 0.01) # 创建时间点 (start,stop,step)
# 调用 ode 对 lorenz 进行求解, 用两个不同的初始值 W1、W2 分别求解
W1 = (0.0, 1.00, 0.0) # 定义初值为 W1
track1 = odeint(lorenz, W1, t)
W2 = (0.0, 1.01, 0.0) # 定义初值为 W2
track2 = odeint(lorenz, W2, t)
# 绘图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot(track1[:,0], track1[:,1], track1[:,2], color='magenta') # 绘制轨
#迹 1
ax.plot(track2[:,0], track2[:,1], track2[:,2], color='deepskyblue') # 绘
#制轨迹 2
ax.set_title("Lorenz attractor by scipy.integrate.odeint")
plt.show()
from scipy.integrate import odeint # 导入 scipy.integrate 模块 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 导数函数, 求 W=[x,y,z] 点的导数 dW/dt def lorenz(W,t): # by youcans p,r,b=10,28,3 x, y, z = W # W=[x,y,z] dx_dt = p*(y-x) # dx/dt = p*(y-x), p: sigma dy_dt = x*(r-z)-y # dy/dt = x*(r-z)-y, r:rho dz_dt = x*y - b*z # dz/dt = x*y - b*z, b;beta return np.array([dx_dt, dy_dt, dz_dt]); t = np.arange(0, 30, 0.01) # 创建时间点 (start,stop,step) # 调用 ode 对 lorenz 进行求解, 用两个不同的初始值 W1、W2 分别求解 W1 = (0.0, 1.00, 0.0) # 定义初值为 W1 track1 = odeint(lorenz, W1, t) W2 = (0.0, 1.01, 0.0) # 定义初值为 W2 track2 = odeint(lorenz, W2, t) # 绘图 fig = plt.figure() ax = Axes3D(fig) ax.plot(track1[:,0], track1[:,1], track1[:,2], color='magenta') # 绘制轨 #迹 1 ax.plot(track2[:,0], track2[:,1], track2[:,2], color='deepskyblue') # 绘 #制轨迹 2 ax.set_title("Lorenz attractor by scipy.integrate.odeint") plt.show()
-
(3)matlab
编辑器:
function dy=vdp1000(t,y)
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
命令窗口:
>> [t,y]=ode15s('vdp1000',[0 3000],[2,0]);
>> plot(t,y(:,1),'-')
-
-
求个关注
数学建模——实验9:微分方程求解
于 2024-05-12 16:29:17 首次发布