0

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

kk75
  • 1

1 Answers1

0

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.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
  • My implementation was a guess,I am meant to implement the implicit scheme as in the two equations given prior to the code,I am not sure how to code this though. If you had not told me this was explicit I would have thought this was correct, though I now realize it is wrong as xnew etc are computed prior to x and y. Could you please show how could I implement an implicit scheme as described in the equations and how could we achieve convergence ? I have only computed explicit schemes and I have tried to look for an example but was not able to find in Python. Do I need 2 for loops? thanks! – kk75 May 04 '20 at 09:14
  • What exactly is $x_{old}$ in your scheme? How is that computed, theoretically? – Lutz Lehmann May 04 '20 at 09:27
  • Sorry for the confusion, I think xold and yold are this last values of $x(t_n+1) $ and $y(t_n+1) $ – kk75 May 04 '20 at 09:52