旨在求解这个耦合微分方程系统。
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例程? 谢谢。
请再研究一下隐式欧拉法或梯形法是如何工作的,对于标量方程和系统。然后仔细想一想你的函数的接口,在这些接口上有一个 x
, y
以及为什么没有 t
在声明中 f
和 g
.
然而,从你的描述来看,你实现的不是隐式方法,而是显式二阶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)
但正如所说,这是一种有固定步数的显式方法,而不是采用隐式方程解策略的隐式方法。