# coupled ode's

Hi everybody,

I am a quite new user of sage, and -up till now- I am really amazed by it power and by its user-friendliness and intuitiveness, I am really happy.

I came up with a problem in solving numerically systems of ode's.

the things are as follows:

I solve a system of ode's (using "desolve_system_rk4") let's say with dependent variables x,y,z and dependent variable t. so, theoretically I end up with 3 functions of time x(t),y(t),z(t). In fact I only have a list of 'points' [t,x(t),y(t),z(t)] for a certain time-grid that I can choose.

now, suppose I have to find many other functions of time (actually a matrix with entries as function of time, let us call it A(t)), whose dynamics (i.e. the coefficients of the ode) depends on x(t) y(t) z(t), how can I give sage the input x(t) y(t) z(t)? I have only a bunch of points.

I tried the obvious way, which is to interpolate the points with "spline", but then the solver of the ode's (for A(t)) doesn't work anymore, it gives me an error of floating precision.

the other way, which actually worked, but to me it seems a bit sloppy, is to solve directly the entire system of x,y,z AND the matrix of functions A(t). This gives the answer, but it seems to me a bit unelegant that I have to recalculate the x,y,z which I already know..can you see the point?

thank you d-

edit retag close merge delete

Sort by » oldest newest most voted

Hi,

I don't know your problem, but sometimes it is possible to introduce the differential equation for A(t) in the set of equations that you pass to the ode solver? Then you don't have to interpolate the solutions.

By the way, you can also try desolve_odeint, which is faster and can deal with compiled code?

Good luck!

Joaquim

more

Indeed 'desolve_odeint' is a much better solver, and the approach Joaquim Puig suggests is the method I usually use for these situations.

( 2015-01-06 14:20:53 +0100 )edit

You can do something like

var('t x y z')
As = [ A.subs( *{ k:v for k,v in zip( (t, x, y, z), pt } ) for pt in list_of_points ]


to get a list of values of $A$, evaluated at each of those points.

You probably need to do that to each entry in the matrix, since subs() doesn't work directly on a whole matrix.

more