ndsolve??

[dabu]

Hi,

I am new in sage. I was wondering about Sage's capability to solve
odes numerically.

I was expecting to find something which is like ndsolve of
Mathematica.
For example it should not only as for the first order equations, nor
that one has to supply jacobians manually. It should also tackle
boundary conditions as equations.

I wonder any such things exist. Without such a construct. Sage would
not be useful (read compete with Mathematica :)) to a large part of
theoretical physics community.

I would be happy to write such a thing in case it does not exist.

thanks,
Pallab

[Jason Grout]

I've always just used the scipy functions to do this sort of thing,
though I'm not an expert in the area, so I'm not sure what the scipy
functionality is missing. See:

http://www.sagenb.org/home/pub/258/ - basic example
http://www.sagenb.org/home/pub/106/ - calculating streamlines in Sage
http://www.sagenb.org/home/pub/879/ - epidemic modeling
http://www.sagenb.org/home/pub/42/ - spring-mass systems

Docs:

http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#ordinary-differential-equations-odeint

http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html#scipy.integrate.odeint

http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode

It would be great to have something like this (or wrap the scipy stuff)
in Sage, so if you want to contribute, go for it! I would be glad to
see it.

Thanks,

Jason

[dabu]

Hi,
Many thanks. I will try to have a look.
best,
Pallab

On Mar 25, 1:05 pm, Jason Grout <jason-s...@creativetrax.com> wrote:
> On 03/25/2010 10:25 AM, dabu wrote:
>
>
>
> > Hi,
>
> > I am new in sage. I was wondering about Sage's capability to solve
> > odes numerically.
>
> > I was expecting to find something which is like ndsolve of
> > Mathematica.
> > For example it should not only as for the first order equations, nor
> > that one has to supply jacobians manually. It should also tackle
> > boundary conditions as equations.
>
> > I wonder any such things exist. Without such a construct. Sage would
> > not be useful (read compete with Mathematica :)) to a large part of
> > theoretical physics community.
>
> > I would be happy to write such a thing in case it does not exist.
>
> I've always just used the scipy functions to do this sort of thing,
> though I'm not an expert in the area, so I'm not sure what the scipy
> functionality is missing.  See:
>
> http://www.sagenb.org/home/pub/258/- basic examplehttp://www.sagenb.org/home/pub/106/- calculating streamlines in Sagehttp://www.sagenb.org/home/pub/879/- epidemic modelinghttp://www.sagenb.org/home/pub/42/- spring-mass systems
>
> Docs:
>

> http://docs.scipy.org/doc/scipy/reference/tutorial/integrate.html#ord...
>
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.o...
>
> http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.o...

[Rajeev]

Hi,
sage has gsl as one of the included packages, which is very good for
numerical solution of differential equations. have a look at examples
on the wikipage -
http://wiki.sagemath.org/interact/diffeq
'Vector Field with Runga-Kutta-Fehlberg' by Schilly is one of my
favorites. i hope it will help.
Best wishes,
Rajeev

On Mar 26, 11:06 am, dabu <pallabb...@gmail.com> wrote:
> Hi,
> Many thanks. I will try to have a look.
> best,
> Pallab
> On Mar 25, 1:05 pm, Jason Grout <jason-s...@creativetrax.com> wrote:
>
> > On 03/25/2010 10:25 AM, dabu wrote:
>
> > > Hi,
>
> > > I am new in sage. I was wondering about Sage's capability to solve
> > > odes numerically.
>
> > > I was expecting to find something which is like ndsolve of
> > > Mathematica.
> > > For example it should not only as for the first order equations, nor
> > > that one has to supply jacobians manually. It should also tackle
> > > boundary conditions as equations.
>
> > > I wonder any such things exist. Without such a construct. Sage would
> > > not be useful (read compete with Mathematica :)) to a large part of
> > > theoretical physics community.
>
> > > I would be happy to write such a thing in case it does not exist.
>
> > I've always just used the scipy functions to do this sort of thing,
> > though I'm not an expert in the area, so I'm not sure what the scipy
> > functionality is missing.  See:
>

> >http://www.sagenb.org/home/pub/258/-basic examplehttp://www.sagenb.org/home/pub/106/-calculating streamlines in Sagehttp://www.sagenb.org/home/pub/879/-epidemic modelinghttp://www.sagenb.org/home/pub/42/-spring-mass systems

[Thierry Dumont]

Rajeev a �crit :

> Hi,
> sage has gsl as one of the included packages, which is very good for
> numerical solution of differential equations. have a look at examples
> on the wikipage -
> http://wiki.sagemath.org/interact/diffeq
> 'Vector Field with Runga-Kutta-Fehlberg' by Schilly is one of my
> favorites. i hope it will help.
> Best wishes,
> Rajeev
>
>

Ok, but... this is a bit sad, but we do not have the best methods for
integrating odes. This would be nice to have (in scipy ?) the most
modern methods, and this would not be very difficult to implement them,
since these method share a rather common interface with odepack (Gear's
method). Real problems are stiff, that is to say (roughly speaking!)
that in dU/dt=F(U), the jacobian of F has eigenvalues whose real part
are distributed in a large negative interval (say... (-10^7, -1) for the
classical example of the Oregonator).

Having a rather long experience in ODE solving, my conclusion is that
the most universal methods are those of the Radau family (described in
the book of Hairer and Wanner "Solving ordinary differential equations"
(part 2). These methods are the most robust of all, but they are also
the fastest ones, for difficult problems.

I wonder whether this would be interesting to have symplectic
integrators, for the integration of Hamiltonian systems.

t.d.

[pallab]

Hi,
@rajeev
Thanks for the gsl link. I am aware of it. gsl is interface is really
horrible :). We need something way smarter.
@dumont
We can surely take into account the stiff solvers. I am interested to
write such codes.

On Mar 28, 2:19 am, Rajeev <rajs2...@gmail.com> wrote:
> Hi,
> sage has gsl as one of the included packages, which is very good for
> numerical solution of differential equations. have a look at examples

> on the wikipage -http://wiki.sagemath.org/interact/diffeq

[Tobias Katz]

When I use substitute_function to declare a function within an existing
equation, I get the following message:

DeprecationWarning: Substitution using function-call syntax and unnamed
arguments is deprecated and will be removed from a future release of
Sage; you can use named arguments instead, like EXPR(x=..., y=...)

Also, when I substitute a function
f(t) with x, x is not interpreted as a constant but as t which is not
what I want.

sage: var('x,t,t1,t2')
(x, t, t1, t2)

sage: function('o',t)
o(t)

sage: integrate(o(t),t,t2,t1)
integrate(o(t), t, t2, t1)

sage: integrate(o(t),t,t2,t1).substitute_function(o,x)
1/2*t1^2 - 1/2*t2^2

(here, I expected something like x*(t1-t2))

What can I do to work around the Warning as well as make sage to do what
I want?

----------------------------------------------------------------------
| Sage Version 4.3.3, Release Date: 2010-02-21 |
| Type notebook() for the GUI, and license() for information. |
----------------------------------------------------------------------
On Ubuntu 9.10

Thank you very much,
Tobi

[pallab]

integrate(f(t),t,t2,t1).substitute({f(t):x})

will do.

[pallab]

Hi,
I am trying to write a little code to create a smart 'ndsolve' routine
like that in mathematica. I will try to interface with gsl ode_solver.
I am less than week old python newbie :), so more experienced people
please help me out...

Let's take the following as an example, defining our function and
equation. Just now I am concentrating on a single equation. After
fixing this simple case, I would try to generalize that later.

f,g,h=function('f',x),function('g',x),function('h',x)
x,z=var('x,z')
eqn2=diff(f(x),x,x)+10*diff(f(x),x)*(f(x)**2-1)+f(x)
bc=[f(0)+2*diff(f,x)(0)==2,diff(f,x)(0)+f(0)==1]

The module to calculate the order of ode (working):

def ode_order2(eqn,f,x,mdeg=10):
L=[]
var('_p')
function('_g')
_g=f
for i in range(1,mdeg):
_g=diff(_g,x)
eqn1=eqn
eqn=eqn.subs({_g:_p})
if (eqn!=eqn1):
L.append(i)
return(max(L))

To transfrom the boundary condition to a array suitable for
ode_solver(), working.

def trans_bc(eqn,bc,f,x,t,mdeg=10):
degree=ode_order2(eqn,f,x,mdeg)
dictarr=[]
p=[var('p'+str(i)) for i in range(degree)]
dictarr.append((f(t),p0))
for i in range(1,degree):
dictarr.append((diff(f(x),*[x for j in range(i)])(t),p[i]))
p_dict=dict(dictarr)
return([j.rhs() for j in solve([i.subs(p_dict) for i in bc],*p)
[0]])

The module to transform the equations (working, but not with
ode_solve()),

def trans_eqn(eqn,f,x,_y,mdeg=10):
dictarr=[]
dictarr.append((f(x),_y[0]))
degree=ode_order2(eqn,f,x,mdeg)
for i in range(degree,0,-1):
dictarr.append((diff(f,*[x for j in range(i)]),_y[i]))
d_dict = dict(dictarr)
yarr=[_y[j] for j in range(degree-1,0,-1)]
sol=solve(eqn.subs(d_dict),_y[degree])[0].rhs()
yarr.append(sol)
return(yarr)

As far as my knowledge goes sage\python do not permit undefined
arrays. So only way to do things is with some named construct,

q=[var('q'+str(i)) for i in range(10)]
trans_eqn(eqn2,f,x,q)

gives,
[q1, -10*(q0^2 - 1)*q1 - q0]

what it should.

But the following code,

def f_1(t,y):
return trans_eqn(eqn2,f(x),x,y)

T=ode_solver()
T.function=f_1
T.y_0=trans_bc(eqn2,bc,f(x),x,0)
T.error_rel=1e-4
T.ode_solve(t_span=[0,1],num_points=100)

is not working.

Traceback (most recent call last): T.y_0=trans_bc(eqn2,bc,f(x),x,0)
File "", line 1, in <module>

File "/tmp/tmpKxA_KC/___code___.py", line 11, in <module>

T.ode_solve(t_span=[_sage_const_0 ,_sage_const_1 ],num_points=_sage_const_100 )
File "", line 1, in <module>

File "ode.pyx", line 531, in sage.gsl.ode.ode_solver.ode_solve (sage/
gsl/ode.c:4558)
ValueError: error solving

thanks,
Pallab