CODEWITHSUNDEEP
pythonmathnumerical-methodsodedifferential-equations
kk75
kk75

Reputation: 1

Unexpected result when solving a system of coupled ode's implicitely

Aiming to solve this system of coupled differential equations:

$ frac{dx}{dt} = -y $

$\frac{dy}{dt} = x $

following the below implicit evolution scheme:

$$ 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})) $$

My code is as follows 

# 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)

I was expecting the output (phase plot) to be a full circle as reported in the literature though I got a spiral which was not expected (I would have embedded the picture though my low reputation did not allowed it, please run to see the output).

Does anyone know how could I fix my Imp routine ? thanks

Upvotes: 0

Views: 61

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

Please study again how the implicit Euler or trapezoidal method works, for scalar equations and for systems. Then think hard about the interfaces of your function, where there is an x, y and why there is no t in the declaration of f and g.

However, from what you describe you are not implementing an implicit method but the explicit 2nd order Heun method. In an implicit method you would solve the implicit equation until sufficient "numerical" convergence is reached.

The explicit Heun loop could look like

    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)

But as said, this is an explicit method with a fixed number of explicit steps, not an implicit method using an implicit equation solving strategy.

Upvotes: 0

Related Questions