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
On Mar 25, 8:25 am, dabu <[email protected]> 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 likendsolveof
> 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
--
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
To unsubscribe, reply using "remove me" as the subject.
