{"action":"created","issue":{"url":"https://api.github.com/repos/brian-team/brian2/issues/626","labels_url":"https://api.github.com/repos/brian-team/brian2/issues/626/labels{/name}","comments_url":"https://api.github.com/repos/brian-team/brian2/issues/626/comments","events_url":"https://api.github.com/repos/brian-team/brian2/issues/626/events","html_url":"https://github.com/brian-team/brian2/issues/626","id":126326003,"number":626,"title":"simple dynamical systems weirdness","user":{"login":"epaxon","id":6740515,"avatar_url":"https://avatars.githubusercontent.com/u/6740515?v=3","gravatar_id":"","url":"https://api.github.com/users/epaxon","html_url":"https://github.com/epaxon","followers_url":"https://api.github.com/users/epaxon/followers","following_url":"https://api.github.com/users/epaxon/following{/other_user}","gists_url":"https://api.github.com/users/epaxon/gists{/gist_id}","starred_url":"https://api.github.com/users/epaxon/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/epaxon/subscriptions","organizations_url":"https://api.github.com/users/epaxon/orgs","repos_url":"https://api.github.com/users/epaxon/repos","events_url":"https://api.github.com/users/epaxon/events{/privacy}","received_events_url":"https://api.github.com/users/epaxon/received_events","type":"User","site_admin":false},"labels":[{"url":"https://api.github.com/repos/brian-team/brian2/labels/component:%20state%20updates","name":"component: state updates","color":"0b02e1"}],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":1,"created_at":"2016-01-13T02:29:15Z","updated_at":"2016-01-13T13:08:43Z","closed_at":null,"body":"Hi,\r\n\r\nI was trying to make some examples of simple 2D linear dynamical systems using brian2. Found some strange behavior/bugs. Here is the code:\r\n\r\n```python\r\nfrom brian2 import *\r\n\r\n# dynamic parameters\r\na = 0.0 / second\r\nb = 1.0 / second\r\nc = -0.5 / second\r\nd = -0.1 / second\r\n\r\neqs_2D = Equations(\"\"\"\r\n    dx/dt = a * x + b * y : 1\r\n    dy/dt = c * x + d * y : 1\r\n\"\"\")\r\n\r\ndefaultclock.dt = 0.1 * ms\r\nrun_time = 40 * second\r\n\r\n# This is just to run several simulation in parallel\r\nn_runs = 10\r\nsystem = NeuronGroup(n_runs, model=eqs_2D, threshold='x>100', reset='x=0; y=0', refractory=0*ms)\r\n\r\n# set the initial conditions\r\nsystem.x = np.random.randn(n_runs)\r\nsystem.y = np.random.randn(n_runs)\r\n\r\nsystem_monitor = StateMonitor(system, True, record=True)\r\n\r\nsimple2D = Network()\r\n\r\nsimple2D.add(system)\r\nsimple2D.add(system_monitor)\r\nsimple2D.store('init')\r\nsimple2D.restore('init')\r\nsimple2D.run(run_time)\r\n```\r\n\r\nSo, this works fine so long as I set [a,b,c,d] such that the system does not have any spirals (stable or unstable). When I set [a,b,c,d] to have spirals (i.e. like it is set right now, which is just damped harmonic oscillator), then the simulation does not work -- every x and y value after the initial condition in the system_monitor is nan. Saddle points and stable stars run fine.\r\n\r\nI've tried a ton of stuff to fix this (changing the dt, changing the parameters), but whenever I get a dynamical system that has a spiral it returns only nan values.\r\n\r\nI finally did something completely trivial to the equations:\r\n\r\n```python\r\neqs_2D = Equations(\"\"\"\r\n    dx/dt = a * x + b * y : 1\r\n    dy/dt = c * x + d * y : 1\r\n    dz/dt = 0 / ms : 1\r\n\"\"\")\r\n```\r\n\r\nand then it works perfectly!\r\n\r\n```python\r\nfigure('x and y vs t')\r\nsubplot(211)\r\nplot(system_monitor.t / second, system_monitor.x.T)\r\nsubplot(212)\r\nplot(system_monitor.t / second, system_monitor.y.T)\r\n\r\nfigure('x vs y')\r\nplot(system_monitor.x.T, system_monitor.y.T)\r\n```\r\n\r\n![image](https://cloud.githubusercontent.com/assets/6740515/12283227/7b613256-b959-11e5-8ef0-a66b8efa02b7.png)\r\n![image](https://cloud.githubusercontent.com/assets/6740515/12283231/813978c8-b959-11e5-9a42-50592b85ede5.png)"},"comment":{"url":"https://api.github.com/repos/brian-team/brian2/issues/comments/171285232","html_url":"https://github.com/brian-team/brian2/issues/626#issuecomment-171285232","issue_url":"https://api.github.com/repos/brian-team/brian2/issues/626","id":171285232,"user":{"login":"mstimberg","id":1381982,"avatar_url":"https://avatars.githubusercontent.com/u/1381982?v=3","gravatar_id":"","url":"https://api.github.com/users/mstimberg","html_url":"https://github.com/mstimberg","followers_url":"https://api.github.com/users/mstimberg/followers","following_url":"https://api.github.com/users/mstimberg/following{/other_user}","gists_url":"https://api.github.com/users/mstimberg/gists{/gist_id}","starred_url":"https://api.github.com/users/mstimberg/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/mstimberg/subscriptions","organizations_url":"https://api.github.com/users/mstimberg/orgs","repos_url":"https://api.github.com/users/mstimberg/repos","events_url":"https://api.github.com/users/mstimberg/events{/privacy}","received_events_url":"https://api.github.com/users/mstimberg/received_events","type":"User","site_admin":false},"created_at":"2016-01-13T13:08:43Z","updated_at":"2016-01-13T13:08:43Z","body":"Hi, thanks for reporting this. I am not sure this strictly qualifies as a bug in our numerical integration, but we should certainly try to make situations like this clearer and less surprising (it is also somewhat related to #554). By default, if you don't specify any numerical integration method for your `NeuronGroup`, Brian tries to solve your equations exactly (over a single time step). This is what it did here, it was able to symbolically solve the equations with `sympy`. Without introducing new constants, the solution is quite complex, though: https://www.wolframalpha.com/input/?i={x%28t%29%27+%3D+a*x%28t%29+%2B+b*y%28t%29%2C+y%28t%29%27+%3D+c*x%28t%29+%2B+d*y%28t%29}\r\n\r\nThe problem with the solution is that depending on the values of the constants, the radicands can become negative -- we don't support complex arithmetics so this will result in not-a-number. Now, why does this change when you introduce `dz/dt = 0 / ms : 1`? In this case, we are no longer able to solve the linear system (which *is* a bug, I guess) so we will use a simple forward Euler integration. \r\n\r\nLong story short: for this kind of system I would not rely on Brian's ability to solve the system explicitly, instead use a numerical solver (specify `euler`, `rk2`, or `rk4` as the `method` argument of `NeuronGroup`). Since your system is linear, `exponential_euler` would also be a good choice but the solution contains divisions by `a` and `d`, so this won't work when these constants are zero. To avoid such problems, it can be a good idea to include the constants directly as values in the equations, instead of as symbols (for now the integration algorithm always leaves constants as symbols in the code, we might change this in the future though: #480). The easiest way to do this is to write your equations as `eqs_2D = Equations(\"\"\"...\"\"\", a=a, b=b, c=c, d=d)`. Writing that, I just found out that this won't work with `exponential_euler` currently due to a misguided, too conservative, check whether the algorithm can integrate the equations..."}}