Неожиданный результат при имплицитном решении системы связанных од - PullRequest
Новая
       137

Неожиданный результат при имплицитном решении системы связанных од

0 голосов
/ 03 мая 2020

Стремление решить эту систему связанных дифференциальных уравнений:

$ 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? спасибо

1 Ответ

0 голосов
/ 04 мая 2020

Пожалуйста, изучите еще раз, как работает неявный метод Эйлера или трапеции, для скалярных уравнений и для систем. Затем подумайте об интерфейсах вашей функции, где есть x, y и почему в объявлении f и g.

нет t. Однако из что вы описываете, вы реализуете неявный метод, а явный метод Heun 2-го порядка. В неявном методе вы решали бы неявное уравнение до тех пор, пока не будет достигнута достаточная «числовая» сходимость.

Явное Heun l oop может выглядеть как 

    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)

Но, как сказано, это явный метод с фиксированным числом явных шагов, а не неявный метод, использующий стратегию неявного решения уравнений.

...