help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

simulating dynamical systems with arbitrary external inputs using lsode


From: John W. Eaton
Subject: simulating dynamical systems with arbitrary external inputs using lsode
Date: Fri, 22 Jun 2007 15:27:28 -0400

On 22-Jun-2007, Scott Kuntze wrote:

| Suppose I have a simple dynamical system
| 
|     function xdot = f(x,t)
|         xdot = 0 ;
|         xdot = -x + du(t)
|     endfunction
| 
| that I wish to simulate with
| 
|     lsode('f',0,t).
| 
| Suppose also that the time vector is
| 
|     t = [0, 0.1, 0.2, 0.3, 0.4]
| 
| and that I have an arbitrary vector, say
| 
|     waveform = [2, 3, 5, 6, 2],
| 
| that I captured from an experiment; the elements in waveform correspond 
| to the
| times in t.
| 
| 
| How can I simulate 'f' with 'waveform' as the input du(t)?  I don't know 
| how to
| get t in the function f to pull out the corresponding element in waveform.
| 
| I believe this is an indexing problem -- I need to link a particular 
| time t1 in
| t to its index, so that the index can then be used to pull out the 
| element in
| waveform at the same index e.g. when t = 0.3, du(0.3) = waveform(4) = 6.
| 
| I tried an ugly solution like the following:
| 
|     global waveform
|     global tt
|     tt = t
|     
|     function result = du(t)
|         global waveform
|         global tt
|         idx = find(tt==t) ;
|         result = waveform(idx) ;    
|     endfunction,
| 
| but lsode stalls -- it doesn't like my "user-defined function," and I 
| suspect
| that the "find" function is the culprit.
| 
| Any suggestions?  I feel this should be an easy thing to do, but I'm stumped
| despite searching the archives and the web at large...

I think something like the following should work:

  function retval = u (t)
    if (t < 5)
      retval = 1;
    elseif (t < 10)
      retval = -1;
    elseif (t < 15)
      retval = 1;
    else
      retval = 0;
    endif
  endfunction

  function xdot = f (x, t)
    xdot = -x + u(t);
  endfunction

  t = 0:0.1:20;
  y = lsode (@f, 0, t);
  plot (t, y);

jwe


reply via email to

[Prev in Thread] Current Thread [Next in Thread]