显式欧拉法求解非线性微分方程组提前终止问题排查

我尝试用显式欧拉法求解如下非线性微分方程组：

x' = f₁(x,y)
y' = f₂(x,y)

已知解对应的曲线需在x-y平面连接初始点(x_initial,y_initial)与(0,1)，但计算得到的曲线在(0.17,0.98)附近提前终止。我尝试调整参数，但仍无法将结果进一步推向(0,1)。起初我认为方程组在终点附近变为刚性ODE，但查阅刚性ODE相关资料后排除了该可能。请问可能的问题是什么？

我的Python代码 #

import math
import numpy as np
import matplotlib.pyplot as plt

q=1

# 定义f1和f2相关的辅助函数
def l(lna,x,y,m,n,xi,yi): 
  return n *m**(-1)*(np.divide((yi**2)*(np.float64(1)-np.power(x,2)-np.power(y,2))*np.exp(3*(lna-lnai)),(y**2)*(1-xi**2-yi**2)))**(-1/n)

def f1 (x,y,l):
  return -3*x + l*np.sqrt(3/2)* y**2+ 3/2 *x*(2*(x**2)+q*(1-x**2-y**2)) 

def f2 (x,y,l): 
  return -l*np.sqrt(3/2) *y*x + 3/2 *y*(2*x**2+q*(1-x**2-y**2))

# 显式欧拉法实现
def e_E(xa,xb,dlna,m,n,xi,yi):          
    N = int(round((lnaf-lnai)/dlna))     
    lna = np.linspace(0, N*dlna, N+1)
    x = np.empty(N+1) 
    y = np.empty(N+1)
    x[0],y[0] = xi,yi
    for i in range(N):
      sd = l(lna[i],x[i],y[i],m,n,xi,yi)
      x[i+1] = x[i] + dlna * f1(x[i],y[i],sd)
      y[i+1] = y[i] + dlna * f2(x[i],y[i],sd)
    return x,y,lna

# 自变量lna的范围
lnai = np.float64(0)
lnaf = np.float64(15)

# 步长
dlna = np.float64(1e-3)

# 初始条件
yi = np.float64(1e-5)
xi = 0

# 计算并绘图
x1,y1,lna1 = e_E(lnai,  lnaf,  dlna,  np.float64(0.1),       np.float64(2), xi, yi)
plt.plot(x1,y1,'b',label = ('x1'))
plt.legend()
plt.grid()
plt.ylabel('y')
plt.xlabel('x')
plt.show()

解曲线对比 #

我的x-y平面解曲线 #

完整/正确解曲线 #

内容的提问来源于stack exchange，提问作者user280016