# What is wrong with this Euler method code in python?

I am testing the ability of my integration method to reproduce results obtained analytically (for various values of the parameters), but the results turn out to be completely wrong and I can't find what's wrong with my code.

As an example, for the parameter values below, $$x_1$$ is supposed to tend towards $$40$$ with time (t). However, when I run this code in PyCharm, $$x_1$$ stops at $$0.1227$$.

The dynamic equation used is

$$\frac{dx_{1}}{dt} = (N - x_1) \mu - \frac{\theta x_1}{k + x_1^2}$$

Any help would be greatly appreciated, thanks !

# Code:

#Parameters:
N = 50
theta= 0.6
µ = 0.6/412.5
k = 1
n = 2
x1 = 0   #x1(0)
listx1 = []  #vector of all x1-values
dt = 0.01    #delta t

from math import pow
import matplotlib.pyplot as plt

t= 3600*24  # time is 1 day
while t>0:
listx1.append(x1)
dx1 = (µ * (N-x1) - (theta * x1) / (k + pow(x1, 2))) * dt
x1 = x1 + dx1
t -= dt

for i in range(500):
print("x1=", listx1[i])

print("x1 final:", listx1[-1])

t = []
for i in range(len(listx1)):
t.extend([i])

plt.plot(t, listx1, 'r')
plt.autoscale(enable=True, axis='both', tight=False)
v = [0, 3600 * 24, 0, 60]
plt.axis(v)
plt.show()

• You are running your problem for a large amount of simulation time. What is the phenomena you're testing? Is it a phenomenon that requires some quantity to be conserve? If so to this last question, Explicit Euler may numerically dissipate the conserved quantity and make you get very incorrect answers for long timeframes. You should consider running this system for a short timeframe, like a few seconds, and see how Explicit Euler does relative to the analytical solution. Also consider posting the analytical solution so others might be able to look into it. – spektr Jan 24 '17 at 18:36

I did not check your code, however, the result you are getting is also verified by scipy.integrate.odeint

from scipy.integrate import odeint

def ode(x, t, N, theta, mu, k):
dxdt = (N - x) * mu - theta * x / (k + x**2)

return dxdt

N = 50
theta= 0.6
mu = 0.6/412.5
k = 1

x0 = 0
t = np.linspace(0, 30, 101)
sol = odeint(ode, x0, t, args=(N, theta, mu, k))


Here is the plot of the solution (note that a steady steady has already been reached for $t=30$)

The steady state value is sol[-1], which equals $0.1227$.

However, you say you expected a value of around $40$. Let's see what happens. At the steady state, it must hold (with $x_{ss}\triangleq \lim_{t\rightarrow \infty} x_1(t)$) $$0 = (N - x_{ss}) \mu - \frac{\theta x_{ss}}{k+x_{ss}^2}$$

Therefore, the steady state solution is the root of 3rd order polynomial, i.e., there should be three possible steady state values for the system. These can be found using sympy as follows.

import sympy as sp

t, N, theta, mu, k = sp.symbols('t, N, theta, mu, k', real = True)
x_ss = sp.symbols('x_ss', real = True)
eq = sp.Eq(0, (N - x_ss) * mu - theta * x_ss / (k + x_ss**2))

sol_sym = sp.solve(eq, x_ss)


sol_sym is a list containing the three roots in a symbolic form. Substituting the numerical values for the parameters it evaluates as follows

[s.subs({N: 50, theta: 0.6, mu: 0.6/412.5, k:1}).n() for s in sol_sym]

> [0.12273605326634 - 2.29e-16*I, 10.2908641151715 + 3.09e-16*I,
> 39.5863998315621 - 7.96e-17*I]


Ignoring the essentially zero imaginary components, it follows that the steady state value can indeed take the value $0.1227$ obtained by the numerical integration, as well as $39.586$ and $10.29$.

Clearly, with the initial condition $x_1(0)=0$, the system reaches the first steady state value. You need to consider a different initial value to reach $40$. For example, running the numerical code with $x_1(0) = 35$, the system reaches this value at around $t=10000$

• Thanks a lot for putting me on the right track. Now, do you know why, when I try running your code on my IDE, I get this error with the sympy.solve function: raise TypeError("cannot determine truth value of\n%s" % self) TypeError: cannot determine truth value of -2*N**3/27 - Nk + N*(kmu + theta)/(3*mu) < 0 ? – Alexeï Jan 25 '17 at 16:10
• @Alexei I don't know why you are getting this error. The code runs fine for me (Sympy version 1.0) – Stelios Jan 25 '17 at 16:24