ASKSAGE: Sage Q&A Forum - Latest question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Fri, 27 Jun 2014 05:55:47 -0500Phase portraits of 2-dimensional systemshttp://ask.sagemath.org/question/9423/phase-portraits-of-2-dimensional-systems/I'm trying to plot solutions to two dimensional ordinary differential equations. I found that Sage makes it easy to plot a vector field and, using `ode_solver()`, I can plot solutions on top of the vector field by specifying an initial condition `y_0` and some range of time to run (`t_span`).
However, this method I'm using seems to be quite ad hoc, as I have to choose the right initial conditions and time span / know a lot about my system in order to plot a nice picture. Let's make this more concrete:
Say I want to draw a nice phase portrait for
$\dot{x} = -y$
$\dot{y} = -x$
First I generate the vector field:
var('x y')
VF=plot_vector_field([-y,-x],[x,-2,2],[y,-2,2])
Then I use `ode_solver()` to plot solutions with initial conditions going around the unit circle:
T = ode_solver()
T.function=lambda t,y:[-y[1],-y[0]]
solutions = []
c = circle((0,0), 1, rgbcolor=(1,0,1))
for k in range(0,8):
T.ode_solve(y_0=[cos(k*pi/4),sin(k*pi,t_span=[0,1],num_points=100)
solutions.append(line([p[1] for p in T.solution]))
This generates the following picture:
![](http://oi46.tinypic.com/27y4v80.jpg)
But if I change run the system for one more unit of time (set `t_span=[0,2]`), the picture gets skewed:
![](http://oi49.tinypic.com/2mfj1i8.jpg)
This makes sense, of course, because the magnitude of the vectors along $y=-x$ get big as you get further away from the origin. Similarly, the trajectory along $y=x$ has trouble getting to the origin because the magnitude of those vectors get very small. Which all makes me think there's a better way to do this. Thoughts?dxvxdSun, 14 Oct 2012 08:02:40 -0500http://ask.sagemath.org/question/9423/Fastest way to call special function (elliptic integral) from cython code for gsl ode_solver()http://ask.sagemath.org/question/11058/fastest-way-to-call-special-function-elliptic-integral-from-cython-code-for-gsl-ode_solver/I am using **sage.gsl.ode.ode_solver** to solve an ode, and overloading the **c_f()** function as detailed [here](http://www.sagemath.org/doc/reference/calculus/sage/gsl/ode.html).
The code is working using an approximate **c_f()** function. For the exact one, I need to use the complete elliptic integrals. These are special functions that are not part of the libc.math library. They seem to be available through a couple packages on sage. I found I could call them by using **sage.all.elliptic_ec(k)** and **sage.all.elliptic_kc(k)**
This is causing the code to slow down by a factor of about 10^3. I know those special functions are going to be expensive, but timing the elliptic_ec() function in a separate cell returns 0 for the time, so these functions aren't *that* slow.
I'm wondering if the problem is mostly just that I'm calling a python function from cython and losing the speed up because of that? Is there a better way to do that? [There are](http://https://www.gnu.org/software/gsl/manual/html_node/Legendre-Form-of-Complete-Elliptic-Integrals.html#Legendre-Form-of-Complete-Elliptic-Integrals) c libraries available with those functions - would it be better to import the files into my project, and call them from the cython code? (not sure how to do that...) Sage seems to have some gsl packages already - does it also have those special functions through gsl? Is there a different package with a faster form of those functions?
Disclaimer: I've played around with sage a bit in the past, but this is the first intensive numerical simulation I've attempted with it, so sage/python/cython are all relatively new to me. I'm more accustomed to Mathematica, Matlab/Octave and c++
Here are (I believe) the relevant bits of code:
%cython
from libc.math cimport pow
cimport sage.gsl.ode
import sage.gsl.ode
include 'gsl.pxi'
cdef class zeeman_acceleration(sage.gsl.ode.ode_system):
... other stuff ...
cdef double coilFieldR(self, double rp, double zp, double coilR): # T
cdef double B0, alpha, beta, Q, k
... some basic arithmetic ...
if rp > coilR/10e4 :
gamma = zp/rp
EC = sage.all.elliptic_ec(k)
EK = sage.all.elliptic_kc(k)
Br = B0*gamma/(self.pi*sqrt(Q))*(EC*(1+alpha**2+beta**2)/(Q-4*alpha)-EK)
return Br
return 0
... other stuff (including c_f() function that calls coilField...
argentum2fFri, 27 Jun 2014 05:55:47 -0500http://ask.sagemath.org/question/11058/ode_solver : unable to convert to floathttp://ask.sagemath.org/question/8714/ode_solver-unable-to-convert-to-float/I am trying to solve a system of 38 first-order differential equations using `ode_solver()`. Among many problems, I have a problem with complex numbers. The functions in question are complex functions (equations include conjugate(function)), and initial contitions are also complex numbers.
The error I get is following:
TypeError: Unable to convert 0.0353553390593274 - 0.0353553390593274*I
to float; use abs() or real_part() as desired
`0.0353553390593274 - 0.0353553390593274*I` is the value of one of the inital conditions. I am rather new to `sage`, so I wonder if `ode_solver()` can deal with complex functions and complex numbers as such?
Thank you!HurinSun, 12 Feb 2012 23:58:01 -0600http://ask.sagemath.org/question/8714/endless loop using ode_solvehttp://ask.sagemath.org/question/7612/endless-loop-using-ode_solve/Im solving a system equation within a for loop, but it never ends
But its ok to evaluate f1(n) outside the loop.... why?
def f1(ll):
def f1(ll):
T=ode_solver()
f = lambda t,y:[y[1],eqa.subs(l=ll,s=y[0],w=y[1],x=t)]
T.function=f
T.y_0=[0,1]
T.algorithm="rk4"
T.t_span=[0,1]
T.h=1
T.ode_solve(num_points=40)
s=T.solution[-1][1][0]-1
return s
for i in range(3) :f1(i)
doesn't works except when f(1) ,otherwise the solver is unable to end. But f(1) or eqa(l=1) is the real solution--- :S
where eqa is
var('l s x w')
eqa=((l^2 + l)*s - 2*w*x)/(x^2 - 1)
in the other hand this works
def g(ll):
Sol=desolve_system_rk4([w,eqa(l=ll)],[s,w],ics=[0,0,1],end_points=[1.0],ivar=x,step=0.01)
return (Sol[len(Sol)-1][1]-1)
for i in range(3) :f1(i)
so.. whats the problem with ode_solver?
ngativSun, 22 Aug 2010 16:26:49 -0500http://ask.sagemath.org/question/7612/how to do a Backward integration using ode_solver?http://ask.sagemath.org/question/7625/how-to-do-a-backward-integration-using-ode_solver/Let's say that i want to integrate this from ti=0 to tf=10
T=ode_solver()
f = lambda t,y:[y[1],y[0]^2-y[1]]
T.function=f
T.y_0=[0,.1]
T.t_span=[0,10]
T.ode_solve(num_points=100)
H= [ [T.solution[i][0] , T.solution[i][1][0] ,T.solution[i][1][1] ] for i in range(len(T.solution)) ]
for i in range(5): H[i]
[0, 0, 0.100000000000000]
[0.10000000000000001, 0.0095163351611087642, 0.090486759458914337]
[0.20000000000000001, 0.018128063040304006, 0.081894953217975602]
[0.30000000000000004, 0.025923510483407335, 0.0741488103498774]
[0.40000000000000002, 0.032983605743555654, 0.06717622415610694]
How i do to integrate backwards, from ti=10 to tf=0?
This doesn't works, because the solution that it gives me is a constant = ti:
T=ode_solver()
f = lambda t,y:[y[1],y[0]^2-y[1]]
T.function=f
T.y_0=[10,.05]
T.t_span=[10,0]
T.ode_solve(num_points=100)
H= [ [T.solution[i][0] , T.solution[i][1][0] ,T.solution[i][1][1] ] for i in range(len(T.solution)) ]
for i in range(5): H[i]
[10, 10, 0.0500000000000000]
[10.0, 10.0, 0.050000000000000003]
[9.9000000000000004, 10.0, 0.050000000000000003]
[9.8000000000000007, 10.0, 0.050000000000000003]
[9.7000000000000011, 10.0, 0.050000000000000003]
Any advice using ode_solver()?ngativMon, 23 Aug 2010 10:25:26 -0500http://ask.sagemath.org/question/7625/