数学建模——实验9：微分方程求解

• 求下列微分方程的数值解，画出解的图形。
• 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),'-')

• 

• 求个关注