What to do with a 'Could not find starting point on curve. Stopping continuation.'?

[ulyssek]

Hi there,

and thanks again for this amazing package. I have trouble using a custom function, and I have this error which I don't know how to handle.

My system is pretty simple, here is a sample of my code:

H_str = '[[  0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,], [  5.,   2.,   1.,   0.,   0.,   0.,   0.,   0.,   0.,   0.,], [ 22.,  17.,  12.,   8.,   4.,   2.,   1.,   0.,   0.,   0.,], [ 41.,  36.,  31.,  26.,  21.,  16.,  12.,   7.,   4.,   2.,], [ 60.,  55.,  50.,  45.,  40.,  36.,  31.,  26.,  21.,  16.,], [ 78.,  74.,  69.,  64.,  59.,  54.,  50.,  45.,  40.,  35.,], [ 96.,  92.,  87.,  82.,  78.,  73.,  68.,  63.,  59.,  54.,], [114., 110., 105., 101.,  96.,  91.,  87.,  82.,  77.,  72.,], [132., 127., 123., 118., 114., 109., 105., 100.,  95.,  91.,], [149., 145., 140., 136., 131., 127., 122., 118., 113., 109.,]]'

params = {
        'tauS' : 0.06,
        'gammaS'  : 1.282,
}

domains = {'x1': [0, 1], 'x2': [0, 1]}

eqs =  {
        'x1': '-x1 * tauS + (1 - x1) * gammaS * H(x1,x2)',
        'x2': '-x2 * tauS + (1 - x2) * gammaS * H(x2,x1)',
}
functions = {}

   
DSargs = dst.args(name='My_model')

DSargs.pars = params
DSargs.fnspecs = functions

DSargs.varspecs = eqs

DSargs.tdomain = [0, 30]
DSargs.xdomain = domains

DSargs.vfcodeinsert_start = """
def H(x,y,n=10):   
      import numpy as np
      ss = np.linspace(0,1,n)
      foo = """+H_str+"""
      eps = 1/(n-1)
      x_mi_index = np.argmax(ss > x)-1
      y_mi_index = np.argmax(ss > y)-1
      x_ma_index = x_mi_index+1
      y_ma_index = y_mi_index+1
      x_mi,x_ma,y_mi,y_ma = x_mi_index*eps,x_ma_index*eps,y_mi_index*eps,y_ma_index*eps
      if x-x_mi<1-(y-y_mi):
            a = (foo[x_ma_index][y_mi_index]-foo[x_mi_index][y_mi_index])/eps
            b = (foo[x_mi_index][y_ma_index]-foo[x_mi_index][y_mi_index])/eps
            c = foo[x_mi_index][y_mi_index]-a*x_mi-b*y_mi
      else:
            b = (foo[x_ma_index][y_ma_index]-foo[x_ma_index][y_mi_index])/eps
            a = (foo[x_ma_index][y_ma_index]-foo[x_mi_index][y_ma_index])/eps
            c = foo[x_ma_index][y_ma_index]-a*x_ma-b*y_ma
      return x*a+b*y+c
    
"""
DSargs.ignorespecial = ['H']


dmModel = dst.Vode_ODEsystem(DSargs)
fp_coord = pp.find_fixedpoints(dmModel, n=4, eps=1e-8)
nulls_x, nulls_y = pp.find_nullclines(dmModel, *list(eqs.keys()), n=3, eps=1e-8,
                                          max_step=0.01, fps=fp_coord)

args1 = (nulls_x[:, 0], nulls_x[:, 1], 'b')
kwargs1 = {'linewidth':2, 'label':'s2 nullcline'}
args2 = (nulls_y[:, 0], nulls_y[:, 1], 'g')
kwargs2 = {'linewidth':2, 'label':'s1 nullcline'}
    
plt.plot(*args1,**kwargs1)
plt.plot(*args2,**kwargs2)


        
plt.xlim(0, 1)
        
plt.ylim(0, 1)
plt.gca().set_aspect('equal', 'box')
        
plt.ylabel('s1')
plt.legend()

The function H is just a linearization of a 10*10 numpy matrix. Computing one point of H takes a few minutes, so I did it for 100 points and stored it in H_str. Now, I used the fact that H is continuous to have an approximation for every value of (x,y).

I know that the nullclines exists, but Pydstool seems to fail to find them.

Any idea of what I could do?

Thanks a lot.

[robclewley]

Hi. You now have dmModel.Rhs as a callable python function. I think you should first convince yourself that your implementation of H has the desirable properties (continuity, smoothness, well-conditioned, etc.) that you expect using some manual investigation. I'd plot a bunch of things related to the whole space and also locally to the FPs and perform some numerical analysis. It's likely that you're making an assumption about H that is false in this implementation, or there is somehow an actual mistake in the coding. These are the usual reasons that continuation is failing.