LSODE Question - External Inputs?

classic Classic list List threaded Threaded
20 messages Options
Reply | Threaded
Open this post in threaded view
|

LSODE Question - External Inputs?

damian.harty
O Learned Denizens,

I'm trying to use Octave to reproduce calculations I've done elsewhere, with a global objective of moving away from Excel and Simulink over time.

The particular problem I'm interested in for now is vehicle ride modelling - my US readers will note that my spelling is British, as am I - which I have previously performed with some success in both Excel and Simulink. If I start with a nice simple system with one degree-of-freedom and linear characteristics, I can write

function xdot = f (x, t)
  m = 100; % kg
  k = 10000; % N/m
  c = 200; % Ns/m
  xdot(1) = x(2);
  xdot(2) =  -k/m*x(1)-c/m*x(2);
endfunction

...and then use lsode to integrate it:

endtime=4.0;
srate=100;
x0 = [0.1; 0];
t = linspace (0, endtime, srate*endtime)';
x = lsode ("xdot", x0, t);

If I use non-zero initial conditions I get the faimliar damped response - this would be a pretty poor car with these parameters but that's neither here nor there.

So far, so good.

Now, what I want to is to introduce an excitation function at the base. I have some function u(t) that is known in advance. I can presume that its derivative, udot(t) will also be known in advance (or at least calculated).

But what I can't really see is how to pass these two vectors into the function "xdot". (Actually it ought to be clear I can't work out how to pass anything into xdot without lsode having a sulk. I want to declare m, k and c externally and have them picked up inside xdot but can't work out how to do that, either)

If I write my equations out on paper they become

xdot_1(t) = x_2(t)
xdot_2(t) = -k/m*( x_1(t) - u(t) ) - c/m*( x_2(t) - udot(t) )

but I'm clearly missing something in terms of how to formulate it inside the xdot function in Octave to subtract the right value of u at any moment in time as xdot contains no sense of time as far as I can see.

If anyone has approached this type of problem in Octave and can give me some pointers I would be most grateful - I'm sure it can be done but am equally sure I am just missing something simple. Should I be using DASSL or ode45, for example?

Damian Harty
Senior Research Fellow - Vehicle & System Dynamics
Coventry University
UK
Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Juan Pablo Carbajal
On Thu, Feb 23, 2012 at 4:42 PM, damian.harty
<[hidden email]> wrote:

> O Learned Denizens,
>
> I'm trying to use Octave to reproduce calculations I've done elsewhere, with
> a global objective of moving away from Excel and Simulink over time.
>
> The particular problem I'm interested in for now is vehicle ride modelling -
> my US readers will note that my spelling is British, as am I - which I have
> previously performed with some success in both Excel and Simulink. If I
> start with a nice simple system with one degree-of-freedom and linear
> characteristics, I can write
>
> function xdot = f (x, t)
>  m = 100; % kg
>  k = 10000; % N/m
>  c = 200; % Ns/m
>  xdot(1) = x(2);
>  xdot(2) =  -k/m*x(1)-c/m*x(2);
> endfunction
>
> ...and then use lsode to integrate it:
>
> endtime=4.0;
> srate=100;
> x0 = [0.1; 0];
> t = linspace (0, endtime, srate*endtime)';
> x = lsode ("xdot", x0, t);
>
> If I use non-zero initial conditions I get the faimliar damped response -
> this would be a pretty poor car with these parameters but that's neither
> here nor there.
>
> So far, so good.
>
> Now, what I want to is to introduce an excitation function at the base. I
> have some function u(t) that is known in advance. I can presume that its
> derivative, udot(t) will also be known in advance (or at least calculated).
>
> But what I can't really see is how to pass these two vectors into the
> function "xdot". (Actually it ought to be clear I can't work out how to pass
> anything into xdot without lsode having a sulk. I want to declare m, k and c
> externally and have them picked up inside xdot but can't work out how to do
> that, either)
>
> If I write my equations out on paper they become
>
> xdot_1(t) = x_2(t)
> xdot_2(t) = -k/m*( x_1(t) - u(t) ) - c/m*( x_2(t) - udot(t) )
>
> but I'm clearly missing something in terms of how to formulate it inside the
> xdot function in Octave to subtract the right value of u at any moment in
> time as xdot contains no sense of time as far as I can see.
>
> If anyone has approached this type of problem in Octave and can give me some
> pointers I would be most grateful - I'm sure it can be done but am equally
> sure I am just missing something simple. Should I be using DASSL or ode45,
> for example?
>
> Damian Harty
> Senior Research Fellow - Vehicle & System Dynamics
> Coventry University
> UK
>
> --
> View this message in context: http://octave.1599824.n4.nabble.com/LSODE-Question-External-Inputs-tp4414228p4414228.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

Hi,
I am not sure that I understand your problem.
Are you trying to actuate your system with an external input u(t)?
If so, then you just add the external input on the corresponding
variable. You can check the function core/nloscillator.m in the
package mechanics to see an example how to do that (there may be other
ways).
If you also need to add the derivative of that function you could also
do it, in principle lsode accepts any right hand-side f(x,t).

If this is not what you meant, I am failing to understand the problem.

Let me know.

--
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Olaf Till
In reply to this post by damian.harty
On Thu, Feb 23, 2012 at 07:42:50AM -0800, damian.harty wrote:
> ...
> But what I can't really see is how to pass these two vectors into the
> function "xdot". (Actually it ought to be clear I can't work out how to pass
> anything into xdot without lsode having a sulk. I want to declare m, k and c
> externally and have them picked up inside xdot but can't work out how to do
> that, either)
> ...

This is an FAQ. Maybe there is already something in the manual for this, if
not, it should be included ...

To pass additional arguments for a passed function (for example with
lsode) you can use anonymous functions, like that:

# a and b are arguments with additional information
a = linspace (0, 1, 1001);
b.x = .1;
b.y = .2;

X = lsode (@ (x, t) fdot (x, t, a, b), x0, t);

If you are not familiar with anonymous functions, you should refer to the
manual first.

In the example, your user function would have to be something like

function ret = fdot (X, T, A, B)
...

A and B will be set to a and b when the function is called by lsode.

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

signature.asc (853 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

RE: LSODE Question - External Inputs?

Allen.Windhorn-2
In reply to this post by damian.harty
Damian,

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of damian.harty

> I'm trying to use Octave to reproduce calculations I've done
> elsewhere, with a global objective of moving away from Excel and
> Simulink over time.
> ...
> Now, what I want to is to introduce an excitation function at the
> base. I have some function u(t) that is known in advance. I can
> presume that its derivative, udot(t) will also be known in advance
> (or at least calculated).
>
> But what I can't really see is how to pass these two vectors into
> the function "xdot". (Actually it ought to be clear I can't work out
> how to pass anything into xdot without lsode having a sulk. I want
> to declare m, k and c externally and have them picked up inside xdot
> but can't work out how to do that, either)

A trick I have used is to increase the size of the X vector and pass
them in as "stealth" members -- their derivatives would be set to zero
so they wouldn't change.  Another way I have used is global variables
(yuck).  It may also be possible to use a nested function for xdot, so
that it shares the same namespace as the calling function.

I suspect there is a more elegant way, so will be watching the
responses.

> If I write my equations out on paper they become

> xdot_1(t) = x_2(t)
> xdot_2(t) = -k/m*( x_1(t) - u(t) ) - c/m*( x_2(t) - udot(t) )

> but I'm clearly missing something in terms of how to formulate it
> inside the xdot function in Octave to subtract the right value of u
> at any moment in time as xdot contains no sense of time as far as I
> can see.

The calling function is explicitly passing t, so you have access to
that within the xdot function.  If your u(t) is defined by a table,
you will have to define a function that interpolates the table, since
the timesteps are variable.

Regards,
Allen
--
Allen Windhorn, P.E. (MN), CEng  (507) 345-2782
Kato Engineering
P.O. Box 8447, N. Mankato, MN  56002
[hidden email]

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

damian.harty
In reply to this post by Olaf Till
I made some progress passing in the variables for m, k and c:

one_dof.m:

m =   100; % kg
k = 10000; % N/m
c =   200; % Ns/m
endtime=4.0;
srate=100;
x0 = [0.1; 0];
t = linspace (0, endtime, srate*endtime)';
x = lsode (@ (x,t) EqMotion (x,t,m,k,c), x0, t);

EqMotion.m:

function xdot = EqMotion (x, t, m, k, c)
  xdot(1) = x(2);
  xdot(2) = -k/m*x(1)-c/m*x(2);
endfunction

I should politely point out that if the impression is widely held that the manual entry for anonymous functions explains that, then it doesn't - your reply was much clearer than the manual entry.

Emboldened by this, I thought I would try passing in my "u" vector. The u vector is a displacement applied to the base of the system, determined in advance - an "obstacle", if you like. The example is a pretty trivial one but it paves the way for various things now.

This was my attempt:

one_dof.m:

m =   100; % kg
k = 10000; % N/m
c =   200; % Ns/m
endtime=4.0;
srate=100;
x0 = [0; 0];
t = linspace (0, endtime, srate*endtime)';
t_vector=t;
u=t*0;
u(2:10)=0.1;
udot=diff(u);
udot(srate*endtime)=udot(srate*endtime-1);
udot=udot*srate;
x = lsode (@ (x,t) EqMotion (x,t,m,k,c,u,udot,t_vector), x0, t);

EqMotion.m:

function xdot = EqMotion (x, t, m, k, c, u, udot, t_vector)
  u_now=interp1(t_vector,u,t);
  udot_now=interp1(t_vector,udot,t);
  xdot(1) = x(2);
  xdot(2) = k/m*(u_now-x(1))+c/m*(udot_now-x(2));
endfunction

So, right now I see what I expect to see. Whether or not it's the best way to code it I have no idea, in particular making a copy of the time vector to pass in twice seems strange but I can't work out how not to need it for interpolation of the u vector I passed in. If there is a better way to take the derivative on the fly I'd like to hear it, too.

Other than that, I think it's mission accomplished for the moment.

Many thanks for both replies, without which I would have failed dismally!

Damian
Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Juan Pablo Carbajal
On Thu, Feb 23, 2012 at 6:46 PM, damian.harty
<[hidden email]> wrote:

> I made some progress passing in the variables for m, k and c:
>
> one_dof.m:
>
>
>
> EqMotion.m:
>
>
>
> I should politely point out that if the impression is widely held that the
> manual entry for anonymous functions explains that, then it doesn't - your
> reply was much clearer than the manual entry.
>
> Emboldened by this, I thought I would try passing in my "u" vector. The u
> vector is a displacement applied to the base of the system, determined in
> advance - an "obstacle", if you like. The example is a pretty trivial one
> but it paves the way for various things now.
>
> This was my attempt:
>
> one_dof.m:
>
>
>
> EqMotion.m:
>
>
>
> So, right now I see what I expect to see. Whether or not it's the best way
> to code it I have no idea, in particular making a copy of the time vector to
> pass in twice seems strange but I can't work out how not to need it for
> interpolation of the u vector I passed in. If there is a better way to take
> the derivative on the fly I'd like to hear it, too.
>
> Other than that, I think it's mission accomplished for the moment.
>
> Many thanks for both replies, without which I would have failed dismally!
>
> Damian
>
> -----
> Senior Research Fellow - Vehicle & System Dynamics
> Coventry University
> United Kingdom
> --
> View this message in context: http://octave.1599824.n4.nabble.com/LSODE-Question-External-Inputs-tp4414228p4414606.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

You could pass function handles instead of vector.
There is the function data2fun.m, that creates a function handle form
data. I think is well documented, if it is not I will be happy to
improve the help. You can get that function from signal-1.1.2. It is
very simple, it uses interp1 to generate the handle. If you check
inside you will see how to do it.

Once generated you pass the handle instead of the vector and inside
your function you evaluate it at the current time. Check the second
demonstration here
http://octave.sourceforge.net/signal/function/sigmoid_train.html it
uses lsode in this way.

Good luck.


--
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

damian.harty
I don't mind the code I have now unless there is some particular reason it isn't preferred?

However, I do have a new quest (same model) which is to pass back some of the intermediate variables.

I've re-organised the code a bit to give myself a spring preload and a gravity force, this allows me to branch when the force gets to zero and let the system go airborne. It works just fine (as I said I am migrating calculations from another environment so I have a reference solution to which I can compare) but I would like to be able to retrieve time histories from inside the function, say for Fk and Fc:

function [x, ISTATE, MSG] = EqMotion (x, t, m, k, c, u, udot, t_vector)
  u_now=interp1(t_vector,u,t);
  udot_now=interp1(t_vector,udot,t);
  xdot(1) = x(2);
  preload = m*9.81;
  Fk = k*( u_now-x(1) ) + preload;
  if (Fk > 0.)
     Fc = c*( udot_now-x(2) );
    else
     Fk = 0;
     Fc = 0;
  endif
  xdot(2) = Fk/m + Fc/m - 9.81;
endfunction

Any suggestions to retrieving Fk & Fc as a function of t from the main routine (below)?

m =   100; % kg
k = 10000; % N/m
c =   200; % Ns/m
endtime=4.0;
srate=200;
x0 = [0; 0];
t = linspace (0, endtime, srate*endtime)';
t_vector=t;
u=t*0;
u(2:50)=0.3;
udot=diff(u);
udot(srate*endtime)=udot(srate*endtime-1);
udot=udot*srate;
[x, ISTATE, MSG] = lsode (@ (x,t) EqMotion (x,t,m,k,c,u,udot,t_vector), x0, t);

Obviously I could recalculate them using x but I was hoping to retrieve them from the integrator, otherwise I have to update the same calculations in two places, which feels like a bad idea. Unless I use (drum roll) a function file, I suppose...

Damian
Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Olaf Till
On Fri, Feb 24, 2012 at 05:40:47AM -0800, damian.harty wrote:
> I don't mind the code I have now unless there is some particular reason it
> isn't preferred?

It seems that you wanted to send some code with your last e-mail,
but it was not there.

> However, I do have a new quest (same model) which is to pass back some of
> the intermediate variables.

Lsode (and other solvers) will call your user function with unpredictable
values of the state variables and for the independent variable. So it
would be not much use to fiddle with global variables to get intermediate
values back. It would be an idea to expand lsode to keep track of
intermediate variables, but currently it does not do so.

But the problem is solved easily. Make your function 'xdot' return more
than one variable. Only the first returned variable, the derivative, will
be used by lsode. In the second and further returned variables, you can
return intermediates from xdot. After lsode is ready, you will have the
values of the state variables corresponding to the supplied values of the
independent variable. Now feed these values into xdot and you will have
the intermediate variables as its second and further output arguments.

For this to work efficiently, you should program xdot in a way that it
is able to handle vector inputs (with vectorized computations), though
lsode itself will only call it with scalars. If not vectorized, you will
have to call xdot in a loop, which is much slower.

(I did not examine the rest of your mail, hope I havn't missed the point.)

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

signature.asc (853 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

damian.harty
Hmm, I'm becoming baffled with the included code not being there thing. I just pasted it in as plain text last time, and I could see it in my e-mail and from inside Nabble...

Anyway, I think you are understanding what I'm trying to do and I get the non-constant calling interval with lsode.

In principle I thought it wouldn't be any bother to get my function to return more than one variable but what with the anonymous function and ISTAGE and MSG variables I can't really work out what the syntax is to actually achieve it.

Inside the file EqMotion.m I have

function [xdot, ISTATE, MSG] = EqMotion (x,t,m,k,c,u,udot,t_vector)
 ...
 Fx=k*(u_now-x(1))+preload
 ...
endfunction

If I add xdot(3)=Fk into it then lsode protests about a dimensional mismatch.

And inside one_dof.m I have

[x, ISTATE,MSG] = lsode (@ (x,t) EqMotion ((x,t,m,k,c,u,udot,t_vector), x0, t);

So how do I modify the function definition and call to pass back Fk? I want to put it on the left hand side but lsode protests again...

Thanks in advance,

Damian


And
Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Generic Mapping Tools

Ismael Diego Nunez-Riboni
In reply to this post by Olaf Till
Anybody combining Octave (3.4.3) with the generic mapping tools of the
Hawaii University (http://gmt.soest.hawaii.edu/)?

I already googled, I'm aware of the MEX files to read grid files under
Octave, but I wonder if there is an easy way of plotting from Octave
using the GMT as a backend, i.e., not necessarily dreaming of this:

graphics_toolkit gmt

but at least some functions in Octave that pass data to be plotted as
arguments to the GMT functions (not by saving on disk) to plot the data
using GMT (and ideally in a X11 window and not to a PS file).

Thanks a lot for the answer in advance... Cheers, Ismael.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Generic Mapping Tools

bpabbott
Administrator

On Feb 24, 2012, at 11:04 AM, Ismael Núñez-Riboni wrote:

> Anybody combining Octave (3.4.3) with the generic mapping tools of the Hawaii University (http://gmt.soest.hawaii.edu/)?
>
> I already googled, I'm aware of the MEX files to read grid files under Octave, but I wonder if there is an easy way of plotting from Octave using the GMT as a backend, i.e., not necessarily dreaming of this:
>
> graphics_toolkit gmt
>
> but at least some functions in Octave that pass data to be plotted as arguments to the GMT functions (not by saving on disk) to plot the data using GMT (and ideally in a X11 window and not to a PS file).
>
> Thanks a lot for the answer in advance... Cheers, Ismael.

I haven't tried it, but are you looking for something like this ?

        http://code.google.com/p/mirone/

Ben

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: Generic Mapping Tools

Ismael Diego Nunez-Riboni
Thanks for the answer.

On 02/24/2012 05:16 PM, Ben Abbott wrote:

>
> On Feb 24, 2012, at 11:04 AM, Ismael Núñez-Riboni wrote:
>
>> Anybody combining Octave (3.4.3) with the generic mapping tools of the Hawaii University (http://gmt.soest.hawaii.edu/)?
>>
>> I already googled, I'm aware of the MEX files to read grid files under Octave, but I wonder if there is an easy way of plotting from Octave using the GMT as a backend, i.e., not necessarily dreaming of this:
>>
>> graphics_toolkit gmt
>>
>> but at least some functions in Octave that pass data to be plotted as arguments to the GMT functions (not by saving on disk) to plot the data using GMT (and ideally in a X11 window and not to a PS file).
>>
>> Thanks a lot for the answer in advance... Cheers, Ismael.
>
> I haven't tried it, but are you looking for something like this ?
>
> http://code.google.com/p/mirone/
>
> Ben
>

More or less, I guess the GDAL is a library for manipulating geospatial
data (probably similar to the climate data operators) but the package
does not seem to be able to plot... I would like a similar interface but
with the GMTs (which are done to plot)... I guess there is nothing like
that... It's not bad, I only need to combine one or two GMT functions
with Octave, but I wanted to know if there was already worked done in
this direction...
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Juan Pablo Carbajal
In reply to this post by damian.harty
On Fri, Feb 24, 2012 at 3:43 PM, damian.harty
<[hidden email]> wrote:

> Hmm, I'm becoming baffled with the included code not being there thing. I
> just pasted it in as plain text last time, and I could see it in my e-mail
> and from inside Nabble...
>
> Anyway, I think you are understanding what I'm trying to do and I get the
> non-constant calling interval with lsode.
>
> In principle I thought it wouldn't be any bother to get my function to
> return more than one variable but what with the anonymous function and
> ISTAGE and MSG variables I can't really work out what the syntax is to
> actually achieve it.
>
> Inside the file EqMotion.m I have
>
> function [xdot, ISTATE, MSG] = EqMotion (x,t,m,k,c,u,udot,t_vector)
>  ...
>  Fx=k*(u_now-x(1))+preload
>  ...
> endfunction
>
> If I add xdot(3)=Fk into it then lsode protests about a dimensional
> mismatch.
>
> And inside one_dof.m I have
>
> [x, ISTATE,MSG] = lsode (@ (x,t) EqMotion ((x,t,m,k,c,u,udot,t_vector), x0,
> t);
>
> So how do I modify the function definition and call to pass back Fk? I want
> to put it on the left hand side but lsode protests again...
>
> Thanks in advance,
>
> Damian
>
>
> And
>
> -----
> Senior Research Fellow - Vehicle & System Dynamics
> Coventry University
> United Kingdom
> --
> View this message in context: http://octave.1599824.n4.nabble.com/LSODE-Question-External-Inputs-tp4414228p4417389.html
> Sent from the Octave - General mailing list archive at Nabble.com.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave

Hi,

I suggest once more that you take a look at nloscillator form the
mechanics package. You will not loose anything by looking and you may
find solutions that my give other ideas. As always, however Olaf gave
an excellent answer.

Regarding sharing code. You could use the not-yet-official platform to
share code.
http://agora.panocha.org.mx/

--
M. Sc. Juan Pablo Carbajal
-----
PhD Student
University of Zürich
http://ailab.ifi.uzh.ch/carbajal/
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

damian.harty
For real systems a polynomial approach as implemented in nloscillator isn't always helpful - I'm aiming for a fully general piecewise numerical implementation. I've done it very successfully in both Excel and Simulink to guide the Mini Countryman World Rally Car suspension development with Öhlins Racing and I'm trying to reproduce it in Octave in order to do even better in future. I love the Jacobian approach generally, it's great for a large number of things - but it's not what I'm trying to do, and in order to be able to debug things I need to be able to access the intermediate calculations.

Hence I'm looking to be able to interrogate some kind of intermediate variable as well as the state vector on exit.

Rather than coding something new, I'm pretty sure I could do it with the existing Octave functionality if I could just get past the syntax issues I'm having with passing an intermediate variable back out. I'm convinced Olaf has understood what I'm trying to do and has given me the answer, but I'm equally convinced I haven't understood it yet...

Damian



Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Olaf Till
In reply to this post by damian.harty
On Fri, Feb 24, 2012 at 06:43:42AM -0800, damian.harty wrote:

> Hmm, I'm becoming baffled with the included code not being there thing. I
> just pasted it in as plain text last time, and I could see it in my e-mail
> and from inside Nabble...
>
> Anyway, I think you are understanding what I'm trying to do and I get the
> non-constant calling interval with lsode.
>
> In principle I thought it wouldn't be any bother to get my function to
> return more than one variable but what with the anonymous function and
> ISTAGE and MSG variables I can't really work out what the syntax is to
> actually achieve it.
>
> Inside the file EqMotion.m I have
>
> function [xdot, ISTATE, MSG] = EqMotion (x,t,m,k,c,u,udot,t_vector)
>  ...
>  Fx=k*(u_now-x(1))+preload
>  ...
> endfunction
>
> If I add xdot(3)=Fk into it then lsode protests about a dimensional
> mismatch.
>
> And inside one_dof.m I have
>
> [x, ISTATE,MSG] = lsode (@ (x,t) EqMotion ((x,t,m,k,c,u,udot,t_vector), x0,
> t);
>
> So how do I modify the function definition and call to pass back Fk? I want
> to put it on the left hand side but lsode protests again...
>
> Thanks in advance,
>
> Damian
Forgive me when I do not dive into the specifities of your example but
try to give a general explanation. First point: EqMotion is supposed
(by lsode) to return only 'xdot', not 'ISTATE' or 'MSG'. But though
lsode only uses 'xdot', you can return additional variables, e.g.:

function [xdot, Fx] = EqMotion (x, t, m, k, c, u, udot, t_vector)

...

Fx = ...;
...
xdot(1) = ...;
xdot(2) = ...;

endfunction

[X, ISTATE, MSG] = \
lsode (@ (x, t) EqMotion (x, t, m, k, c, u, udot, t_vector), x0, T);

Then, call

[xdot_not_needed_now, Fx] = EqMotion (X.', T, m, k, c, u, udot, t_vector);

to get the intermediate result Fx. For this to work, EqMotion must
be able to handle matrices of X and vectors of T. Otherwise,
you must use an ugly loop:

Fx = zeros (size (T));
for id = 1 : length (T)
  [xdot_not_needed_now, Fx(id)] = \
    EqMotion (X(id, :).', T(id), m, k, c, u, udot, t_vector);
endfor

Olaf

--
public key id EAFE0591, e.g. on x-hkp://pool.sks-keyservers.net

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave

signature.asc (853 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Generic Mapping Tools

Rudolf Widmer-Schnidrig
In reply to this post by Ismael Diego Nunez-Riboni
On 24.02.12 17:04, Ismael Núñez-Riboni wrote:

> Anybody combining Octave (3.4.3) with the generic mapping tools of the
> Hawaii University (http://gmt.soest.hawaii.edu/)?
>
> I already googled, I'm aware of the MEX files to read grid files under
> Octave, but I wonder if there is an easy way of plotting from Octave
> using the GMT as a backend, i.e., not necessarily dreaming of this:
>
> graphics_toolkit gmt
>
> but at least some functions in Octave that pass data to be plotted as
> arguments to the GMT functions (not by saving on disk) to plot the
> data using GMT (and ideally in a X11 window and not to a PS file).
>
> Thanks a lot for the answer in advance... Cheers, Ismael.
> _______________________________________________
> Help-octave mailing list
> [hidden email]
> https://mailman.cae.wisc.edu/listinfo/help-octave
>
Dear Ismael,

my understanding is that the upcoming release of GMT (GMT5) will go
along way in the direction you look for.

Ask the GMT mail list at gmt.soest.hawaii.edu for details.            -Ruedi


--
+---------------------------------------------------------------+
               Rudolf Widmer-Schnidrig
               Black Forest Observatory (BFO)
               Heubach 206
       D-77709 Wolfach, Germany
               e-mail: [hidden email]
               phone:  +49 7836 2151
               fax:  +49 7836 955240
               48d 19' 48.24" North  8d 19' 23.13" East

+---------------------------------------------------------------+

_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

damian.harty
In reply to this post by Olaf Till
Ugly loop or not, it works just fine, thanks for your help - I'll be using this method from now on until I learn more.

Damian
Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

Jordi Gutiérrez Hermoso-2
In reply to this post by Olaf Till
On 24 February 2012 16:05, Olaf Till <[hidden email]> wrote:

> [xdot_not_needed_now, Fx] = EqMotion (X.', T, m, k, c, u, udot, t_vector);

I'm not sure if you're aware or not, but if you're using Octave 3.4.x
or newer, you can use the syntax

    [~, out] = foo(bar);

to omit output variables. Some functions will detect that you are
omitting an output variable, so won't attempt to compute it.

- Jordi G. H.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

damian.harty
I'm running 3.4.0 on OSX 10.6.8 and the tilde (~) trick worked just fine.

Many thanks for your collective assistance in getting me off the mark.

Damian
Senior Research Fellow - Vehicle & System Dynamics Coventry University United Kingdom
Reply | Threaded
Open this post in threaded view
|

Re: LSODE Question - External Inputs?

StandardOctaveUser
Hello Damian,

I assume the issue is long resolved, but
I was wondering if you are interested in
comparing results for such codes,

I have a bunch spreading between
multi degree of freedom output generators
to Enhanced Frequency Domain Decomposition,
Structural Health Monitoring codes.

I was looking for people interested to
get these and similar codes faster and
more efficient.

But be warned some codes  are
extremely slow... I have learning vectorization
in my agenda at a later date due to that.