隐式求解耦合颂歌系统时的意外结果

Question

旨在解决这个耦合微分方程组：

$frac｛dx｝｛dt｝=-y$

$\frac｛dy｝｛dt｝=x$

遵循以下隐式进化方案：

$$y(t_{n+1}(=y(t_｛n｝(+\frac｛\Delta t｝｛2｝(x_｛old｝(t_{n+1｝(+x(t_

$$x(t_{n+1}(=x(t_｛n｝(-\frac｛\Delta t｝｛2｝(y_｛old｝(t_{n+1｝(+y(t_

我的代码如下

# coupled ODE's to be solved 
def f(x,y):
return -y
def g(x,y):
return x
#implicit evolution scheme 
def Imp(f,g,x0,y0, tend,N):
t = np.linspace(0.0, tend, N+1)
dt = 0.1 
x1 = np.zeros((N+1, ))
y2 = np.zeros((N+1, ))
xold = np.zeros((N+1, ))
yold = np.zeros((N+1, ))
xxold = np.zeros((N+1, ))
yyold = np.zeros((N+1, ))
x1[0] = x0 
y2[0] = y0
for n in range(0,N):
xold = f(t[n+1], x1[n])
xxold = f(t[n+1], x1[n+1] + xold)
yold = g(t[n], y2[n])
yyold = g(t[n], y2[n+1]+yold)

y2[n+1] = y2[n] + (x1[n]+xxold)*0.5*dt
x1[n+1] = x1[n] - (y2[n]+ yyold)*0.5*dt
return t,x1,y2
t, y1,y2 = Imp(f,g,np.sqrt(2),0.0, 100, 1000)
plt.plot(y1,y2) 

我原以为输出(相位图(会像文献中报道的那样是一个完整的圆圈，尽管我得到了一个意想不到的螺旋(尽管我的低声誉不允许，但我会嵌入图片，请运行查看输出(。

有人知道我该如何修复我的Imp程序吗？感谢

Answer 1

请再次研究隐式欧拉或梯形方法如何适用于标量方程和系统。然后仔细考虑函数的接口，其中有x、y，以及为什么在f和g的声明中没有t。

然而，根据您的描述，您并没有实现隐式方法，而是实现显式的二阶Heun方法。在隐式方法中，您将求解隐式方程，直到达到足够的"数值"收敛。

显式的Heun循环可能看起来像

for n in range(0,N):
xold = x[n] + f(x[n],y[n])*dt
yold = y[n] + g(x[n],y[n])*dt
xnew = x[n] + f(xold, yold)*dt
ynew = x[n] + g(xold, yold)*dt
x[n+1] = 0.5*(xold+xnew)
y[n+1] = 0.5*(yold+ynew)

但正如所说，这是一种具有固定显式步骤数的显式方法，而不是使用隐式方程求解策略的隐式方法。