避免绘制 ODE 不同的解 ODEint答案

【问题标题】：Avoiding plotting ODEs divergent solutions ODEint避免绘制 ODE 不同的解 ODEint
【发布时间】：2020-12-17 20:55:16
【问题描述】：

我正在尝试绘制一个相平面，我希望它看起来不错。然而，由于初始条件，方程组的一些解会发散。有没有什么方法可以让我在解决方案出现分歧时按顺序创建一个 try/except 链，它不会绘制它。这是我的代码：

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
import pylab as pl 

def aux_func(x):  
  y = x[0]-x[1]
  if (np.abs(y) <= 1):
    f = y**3 + 0.5*y 
  else: 
    f = 2*y - np.sign(y)
  return f 

def function(x,t): 
  x1_dot = x[1]  
  x2_dot = -x[1] - aux_func(x)  
  return [x1_dot,x2_dot]

ts = np.linspace(0, 20, 300)  
ic_1 = np.linspace(-1,1,10)
ic_2 = np.linspace(-1,1,10) 

for r1 in ic_1:
  for r2 in ic_2:
    x0 = (r1,r2)
    try:
      xs = odeint(function, x0, ts)
      plt.plot(xs[:,0], xs[:,1],"r-",linewidth=.8) 
    except: 
      pass

# Nombre de los ejes, limites,
plt.xlabel("$x_1$", fontsize=12)
plt.ylabel("$x_2$", fontsize=12)
# plt.tick_params(labelsize=10)
# plt.xticks(np.linspace(0,1,11))
# plt.yticks(np.linspace(0,1,11))
plt.xlim(-1, 1)
plt.ylim(-1, 1)
# Grafica el campo vectorial
X1, X2 = np.mgrid[-1:1:20j,-1:1:20j]
u=X2
d= X1-X2

t = np.zeros(np.shape(d))
for i in range(len(d)):
    for j in range(len(d[0])):
        if np.abs(d[i][j]) > 1:
            t[i][j]= 2*d[i][j]-0.5*np.sign(d[i][j])
        else:
            t[i][j] =d[i][j]**3 + 0.5*d[i][j]
v=-X2-t

pl.quiver(X1, X2, u, v, color = 'b',width = .002)
plt.grid()
plt.title('Plano de Fase Punto 1')
#plt.savefig('FasePunto4.png')
plt.show()

代码正在绘制以下内容：

感谢您的帮助。

【问题讨论】：

标签： python numpy visualization ode

【解决方案1】：

这可以通过完全避免错误的分歧来解决，这样就不需要异常处理了。

这是一个不连续的 ODE，它会导致异常的效果，例如滑模。一种快速解决此问题的方法是通过实现一个混合区域来缓和跳跃，在该区域中，矢量场快速但连续地从一个阶段变化到另一个阶段（有关其他通用解决方法，请参阅Unsure about how to use event function in Matlab）。对此的更改可以实现为

def aux_func(x):  
  def softsign(u): return np.tanh(1e4*u)
  y = x[0]-x[1]
  h = 0.5*(1+softsign(y**2-1)
  # h is about zero for |y|<1 and about 1 for |y|>1
  f1 = y**3 + 0.5*y          # for |y|<1
  f2 = 2*y - softsign(y)     # for |y|>1, note the second mollification
  return (1-h)*f1+h*f2

无需对代码进行进一步更改，即可得到情节

请注意pylab 已过时，它的所有功能也可以通过plt=matplotlib.pyplot 访问。

【讨论】：