Python solve delay differential equations conditionally

Question

I am using dde23 of pydelay package to solve a delay differential equation. My question: How to code a equation conditionally? For example the target equation has two options:

when x>1, dx/dt=0.25 * x(t-tau) / (1.0 + pow(x(t-tau),10.0)) -0.1*x
otherwise, dx/dt=0.25 * x

I tried two approaches, but it seems like neither one worked:

1. Approach 1 did not complain, but it the if else statement has not been interpreted.

2. Approach 2 generated the following errors:

   Found executable c:\mingw\bin\g++.exe c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp: In function 'double f(double, double)': c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:734: error: 'x' cannot be used as a function c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:734: error: 'x' cannot be used as a function c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp: In function 'PyObject* compiled_func(PyObject*, PyObject*)': c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:878: warning: comparison between signed and unsigned integer expressions c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:883: warning: comparison between signed and unsigned integer expressions c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:774: warning: unused variable 'Nhistx_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:775: warning: unused variable 'Shistx_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:776: warning: unused variable 'Dhistx_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:785: warning: unused variable 'NVhistx_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:786: warning: unused variable 'SVhistx_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:787: warning: unused variable 'DVhistx_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:796: warning: unused variable 'NThist_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:797: warning: unused variable 'SThist_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:798: warning: unused variable 'DThist_ar' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:817: warning: unused variable 'Ndiscont' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:818: warning: unused variable 'Sdiscont' c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp:819: warning: unused variable 'Ddiscont' Traceback (most recent call last): File "C:\Documents and Settings\thao\Desktop\mackey-glass.py", line 33, in dde.run() File "C:\Python27\lib\site-packages\pydelay_dde23.py", line 1120, in run compiler = 'gcc') File "C:\Python27\lib\site-packages\scipy\weave\inline_tools.py", line 355, in inline **kw) File "C:\Python27\lib\site-packages\scipy\weave\inline_tools.py", line 482, in compile_function verbose=verbose, **kw) File "C:\Python27\lib\site-packages\scipy\weave\ext_tools.py", line 367, in compile verbose = verbose, **kw) File "C:\Python27\lib\site-packages\scipy\weave\build_tools.py", line 272, in build_extension setup(name = module_name, ext_modules = [ext],verbose=verb) File "C:\Python27\lib\site-packages\numpy\distutils\core.py", line 186, in setup return old_setup(**new_attr) File "C:\Python27\lib\distutils\core.py", line 169, in setup raise SystemExit, "error: " + str(msg) distutils.errors.CompileError: error: Command "g++ -O2 -Wall -IC:\Python27\lib\site-packages\scipy\weave -IC:\Python27\lib\site-packages\scipy\weave\scxx -IC:\Python27\lib\site-packages\numpy\core\include -IC:\Python27\include -IC:\Python27\PC -c c:\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.cpp -o c:\docume~1\thao\locals~1\temp\thong\python27_intermediate\compiler_a77d1132635f0379270bcb96a5e542fc\Release\docume~1\thao\locals~1\temp\thong\python27_compiled\sc_f68f7a878bf7b27c6f72c9e771ec4d311.o" failed with exit status 1 [Finished in 5.8s with exit code 1]

Approach 1 (use the if else statement to update eqns which is a python dict):

import numpy as np
import pylab as pl
from pydelay import dde23

eqn_1a='0.25 * x(t-tau) / (1.0 + pow(x(t-tau),10.0)) -0.1*x'
eqn_1b='0.45 * x'

eqns = { 'x' : eqn_1a if 'x>1' else eqn_1b}

dde = dde23(eqns=eqns, params={'tau': 15})
dde.set_sim_params(tfinal=1000, dtmax=1.0, AbsTol=10**-6, RelTol=10**-3)

histfunc = {'x': lambda t: 0.5 } 
dde.hist_from_funcs(histfunc, 51)
dde.run()

sol1 = dde.sample(2, 3, 0.1)
x1 = sol1['x']

Approach 2 (provide c code to the solver):

eqns = { 'x' : 'f(x, tau)'}

# We can define a c function to be used in the equations
mycode = """
    double f(double x, double tau) {
        if (x>1){
            return (0.25 * x(t-tau) / (1.0 + pow(x(t-tau),10.0)) -0.1*x);
        }
        else{
            return (0.45 * x);
        }
     }
    """
dde = dde23(eqns=eqns, params={'tau': 15}, supportcode=mycode)

dde.set_sim_params(tfinal=1000, dtmax=1.0, AbsTol=10**-6, RelTol=10**-3)

histfunc = {'x': lambda t: 0.5 } 
dde.hist_from_funcs(histfunc, 51)
dde.run()

sol1 = dde.sample(1, 300, 1)
x1 = sol1['x']

Answer 1

The second approach works with one minor change. The value of the delayed variable has to be passed as an argument to the function directly.

eqns = { 'x' : 'f(x, x(t - tau))'}

mycode = """
double f(double x, double x_tau) {
    if (x > 1.0){
        return 0.25 * x_tau / (1.0 + pow(x_tau, 10.0)) -0.1*x;
    }
    else{
        return 0.45 * x;
    }
}
"""

Although this does run, note that it cannot give a very accurate solution. This is due to the discontinuity at x=1. More advanced solvers can explicitly handle such discontinuities (see e.g. http://www.radford.edu/~thompson/ffddes/).

If you want to stick with the pydelay solver for convenience or just to get an overview, you can try to set the maximum step size dtmax small enough to reduce the error.