时间戳

首页 > 解决方案 > 隐式求解耦合 ode 系统时的意外结果

python - 隐式求解耦合 ode 系统时的意外结果

问题描述

旨在求解该耦合微分方程组：

$ 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_{n} )) $$

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

我的代码如下

# 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 例程？谢谢

标签： pythonmathnumerical-methodsodedifferential-equations

解决方案

请再次研究隐式欧拉或梯形方法如何工作，用于标量方程和系统。然后仔细想想你的函数的接口，在哪里有一个x，y以及为什么在和t的声明中没有。fg

但是，根据您的描述，您没有实现隐式方法，而是实现了显式的二阶 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)

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

推荐阅读