How do I make sage to solve a ODE that includes a time-dependent function?

[Jose Guzman]

Hi everybody,

I am starting to learn how to solve ODEs analytically/numerically with Sage, which to my understanding, implies using Maxima functions. I followed the docu here:  http://maxima.sourceforge.net/docs/manual/en/maxima_22.html.

Now, i wanted to solve analytically the following equation:
cm*dV/dt =  iM - gL*(v(t) - EL) , where cm, iM, gL and EL are parameters.

This works nicely assuming that these parameters are constant. However, I would like to solve the equation for iM being a  time-dependent function (for example a square/sin pulse). Something like this:

cm*dV/dt =  iM(t) - gL*(v(t) - EL), where iM(t) is the time-dependent function

# define independent/dependent variables
sage: t=var('t')
sage: v=function('v',t)

# equation parameters
sage: iM, gL, EL, cm = var('iM, gL, EL, cm')

# define equation
sage: eq = diff(v,t) - 1/cm*(iM-gL*(v-EL) ) == 0

# solve it analytically
sage: desolve(de=eq, dvar=v, ivar=t, ics=[0,-65])

Then Sage returns this:-> (EL*gL*e^(gL*t/cm) + (e^(gL*t/cm) - 1)*iM - (EL +65)*gL)*e^(-gL*t/cm)/gL

which is correct.

If I now define a function of the form:

def pulse(tonset, tdur, amp):
    """ returns a square pulse as a function of t, f(t)
    the pulse is defined as follows:
    tonset -- start of pulse
    tdur   -- duration of pulse
    amp    -- amplitude of pulse
    """
    t=var('t') # do i need this?
   
    f(t)= amp*(sign(t-tonset)/2-sign(t-tonset-tdur)/2)
    return f

I tried to define iM as a function of t
sage: mypulse = pulse(tonset =5, tdur = 10, amp =2 )
 
and redefined the equation (not sure if I wrote correctly
sage:  eq = diff(v,t) - 1/cm*(mypulse - gL*(v-EL) ) == 0
sage:  eq = diff(v,t) - 1/cm*(diff(mypulse,t) - gL*(v-EL) ) == 0 # this is probably wrong

In either case, Sage returns the same error:
TypeError: unable to make sense of Maxima expression 'v(t)=-e^-(t*gL/cm)*(at(integrate((EL*gL+signum(t-5)-signum(t-13))*e^(t*gL/cm),t),[t=0,v(t)=-65])-integrate((EL*gL+signum(t-5)-signum(t-13))*e^(t*gL/cm),t)+65*cm)/cm' in Sage

What I am missing here. I got a similar problem with desolve_laplace. Should I use another ODE solver?

I got it running with Mathematica (using DSolve), but it would be great to have it with Sage!

Thanks a lot in advance

Jose.

[Nils Bruin]

On Aug 5, 8:41 am, Jose Guzman <sjm.guz...@googlemail.com> wrote:
> In either case, Sage returns the same error:
> TypeError: unable to make sense of Maxima expression
> 'v(t)=-e^-(t*gL/cm)*(at(integrate((EL*gL+signum(t-5)-signum(t-13))*e^(t*gL/cm),t),[t=0,v(t)=-65])-integrate((EL*gL+signum(t-5)-signum(t-13))*e^(t*gL/cm),t)+65*cm)/cm'
> in Sage

That seems a bug to me. It is probably because maxima's "at" function
has not been translated yet. This should not be hard to fix. In Sage,
similar functionality is available in the form of "substitute" and
"call". Here is a more direct way of arriving at a similar error
(requires 4.7.1 for interfacing with the LISP below, but that should
not be essential for exposing the problem):

sage: M=sage.calculus.calculus.maxima
sage: E=sage.libs.ecl.ecl_eval
sage: expr=E("#$at(derivative(f(x),x),[x=0])$")
sage: M(expr)
?%at('diff(f(x),x,1),[x=0])
sage: M(expr).sage()
[...]
TypeError: unable to make sense of Maxima expression 'at(diff(f(x),x,
1),[x=0])' in Sage

You can help sage by filing a bug report.

[kcrisman]

On Aug 5, 12:07 pm, Nils Bruin <nbr...@sfu.ca> wrote:
> On Aug 5, 8:41 am, Jose Guzman <sjm.guz...@googlemail.com> wrote:
>
> > In either case, Sage returns the same error:
> > TypeError: unable to make sense of Maxima expression
> > 'v(t)=-e^-(t*gL/cm)*(at(integrate((EL*gL+signum(t-5)-signum(t-13))*e^(t*gL/ cm),t),[t=0,v(t)=-65])-integrate((EL*gL+signum(t-5)-signum(t-13))*e^(t*gL/c m),t)+65*cm)/cm'
> > in Sage
>
> That seems a bug to me. It is probably because maxima's "at" function

Hmm. Maxima's at function IS translated, as you can see in that the
final error has 'at' in it. In calculus.py:

delayed_functions = maxima_qp.findall(s)
if len(delayed_functions) > 0:
for X in delayed_functions:
if X == '?%at': # we will replace Maxima's "at" with
symbolic evaluation, not an SFunction
pass
else:
syms[X[2:]] = function_factory(X[2:])
s = s.replace("?%","")

Probably the problem is that this was too naive. The thing you did
with M(expr).sage() bypasses this, because this 'at' is not
universally imported into the global namespace. Maybe it should be?

- kcrisman

[Jose Guzman]

I would be willing to 1) submit a bug report,  and even 2) try to solve it if somebody would assist me with whole process of submitting a track (maybe in #sagemath channel on freenode. That would be a fantastic opportunity to collaborate with the Sage development. I would be very happy to have the option of integrating this kind of functions in Sage.

[kcrisman]

>
> >> You can help sage by filing a bug report.
>
> I would be willing to 1) submit a bug report,  and even 2) try to solve
> it if somebody would assist me with whole process of submitting a track
> (maybe in #sagemath channel on freenode. That would be a fantastic
> opportunity to collaborate with the Sage development. I would be very
> happy to have the option of integrating this kind of functions in Sage.

Great! We always welcome new help of any kind :)

You should definitely go to http://trac.sagemath.org/sage_trac/ and
follow instructions to get an account there, then fill out a new
defect request ticket with the information you have in your original
post (and/or some of the information in my and Nils' posts, as
relevant). A few different examples (such as the desolve_laplace one)
would be helpful too.

As I point out, I think this is indeed there (see
http://trac.sagemath.org/sage_trac/ticket/385), it's just somehow not
connecting with the output here for some reason. You could start
finding it by following the code for desolve - you could put in some
print statements after the 'soln' is found to see what it looks
like.

I think that here the problem is probably that taylor and laplace
coerce into SR, which forces the 'at' substitution, while the
desolvers all use soln.sage(), which bypasses this. We need to unify
our Maxima conversions (and get a lot of them out of calculus.py,
where they don't really belong) but it's been a back-burner issue
since it's usually easier to do little fixes here and there.

- kcrisman

[Jose Guzman]

On 05/08/11 22:00, kcrisman wrote:
>
> Great! We always welcome new help of any kind :)
>
>
>
>

Here my (first) ticket submission:
http://trac.sagemath.org/sage_trac/ticket/11653

hope its OK

I will try to have a look carefully to that code, to see if I could
provide any help!

Thanks!

Jose