Re: lsode and symbolic function handles

[Vikram Garg]

Thanks Colin and Olaf. I will try the function handle approach.

I have another question related to lsode, is there anyway I can access the solution from an older time step when lsode calls the user defined RHS function ?

Basically, my RHS at a time t(n) depends on values of the solution variable at time t(n-1).

Thanks.

Vikram Garg

vikramvgarg.github.io/

On Sat, Feb 29, 2020 at 4:23 AM Olaf Till <address@hidden> wrote:

On Fri, Feb 28, 2020 at 04:34:41PM -0600, Vikram Garg wrote:
> Hello,
>          I would like to be able to evaluate symbolic function handles to
> construct the forcing vector of an ODE.
>
> For example, in the code below,
>
> ## Construct symbolic vector b
> syms b1 b2;
>
> b(1) = b1;
> b(2) = b2;
>
> ## Handle for symbolic vector b
> bh = function_handle(b);
>
> ## function to be passed to lsode
> function xdot = f(x,t)
>
> xdot = zeros(2,1);
>
> xdot(1) = bh(x(1), x(2))*[1;0];
> xdot(2) = bh(x(1), x(2))*[0;1];
>
> endfunction
>
> ## Initial conditions and time vector for lsode
> x0 = [0;1];
>
> t = linspace (0, 500, 1000);
>
> y = lsode("f",x0, t);
>
> Running this gives an error,
>
>
>
> *error: 'bh' undefinederror: lsode: evaluation of user-supplied function
> failed*
>
> Trying to declare the function handle as global (global bh =
> function_handle(b);) did not help.

This is not specific to lsode. If you want to use a global inside a
function, it must not only be declared global outside the function,
but inside the function, too.

But globals shouldn't be used for this. Give the function handle
('bh' in your case) as an argument to 'f':

function xdot = f(x, t, f_handle)

xdot = zeros(2,1);

xdot(1) = f_handle (x(1), x(2))*[1;0];
xdot(2) = f_handle (x(1), x(2))*[0;1];

endfunction

and wrap 'f' into an anonymous function to pass it to lsode (I can't
say anything, though, to your handling of the symbolic vector since I
don't use Octave as a CAS):

bh = ... some function handle ...;

y = lsode(@ (x, t) f (x, t, bh), x0, t);

Olaf

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