106.sagews

# Declare our variables (x and y)
var('x,y')

# import tools we'll need below
from scipy.integrate import odeint
import numpy

# These are our derivatives; f = dx/dt, g = dy/dt

f(x,y)=-y
g(x,y)=x

# A harder set of functions.
# (there is some sort of problem with this more complicated expression 
# with starting point (0,0) and starting point (0,0.5) using the scipy solver.  
# It's some sort of fortran interface issue, I think, and I think 
# it has to do with unicode, maybe.
#f(x,y) = y+exp(x/10)-1/3*((x-1/2)^2+y^3)*x-x*y^3
#g(x,y) = x^3-x+1/100*exp(y*x^2+x*y^2)-0.7*x

# For performance.  This basically constructs a floating point approximation 
# function that is very fast for repeated evaluation.
ff = fast_float(f,'x','y')
gg = fast_float(g, 'x','y')
fx = fast_float(f.diff(x), 'x','y')
fy = fast_float(f.diff(y), 'x','y')
gx = fast_float(g.diff(x), 'x','y')
gy = fast_float(g.diff(y), 'x','y')

# Initial seed point for streamline
initial_value = (0,1)

# the ending value for the streamline (we go to about 2*pi)
t_end = 6.28

# The number of points in the streamline
num_points = 300

# Here we set up the function, the jacobian, and a sequence of time values for which we will approximate the solution.
# odeint is an ODEPACK (lsoda) based solver

# Basically, think of yp as a coordinate pair (x,y), and we are returning 
# the parametric function output (f(x,y), g(x,y)).  We ignore the time value, since our function output does not depend on the time.
def func(yp, t):
    return [ff(yp[0], yp[1]), gg(yp[0], yp[1])]

# Given an yp=(x,y), return the rows of the jacobian matrix.
def jacobian(yp, t):
    return [[fx(yp[0], yp[1]), fy(yp[0],yp[1])], [gx(yp[0], yp[1]), gy(yp[0], yp[1])]]

# This creates a numpy array of t-values between 0 and t_end with num_points entries.
t=numpy.linspace(start=0, stop=t_end, num=num_points).astype(float)

%time
streamline = odeint(func=func, y0=initial_value, t=t, Dfun = jacobian)
print streamline[0:10] # First 10 entries
show(line(streamline),aspect_ratio=1)

[[ 0. 1. ] [-0.02100179 0.99977944] [-0.04199432 0.99911785] [-0.06296833 0.99801553] [-0.08391457 0.99647295] [-0.10482379 0.99449081] [-0.12568676 0.99206998] [-0.1464943 0.98921151] [-0.16723721 0.98591669] [-0.18790635 0.98218695]]

CPU time: 0.72 s, Wall time: 0.76 s

%time
# Now we plot multiple streamlines using scipy:
p=[] # p will be a list of streamlines
for i in [0..1,step=0.1]:
    p.append(line(odeint(func=func, y0=(0, float(i)), t=t, Dfun = jacobian)))
    print i
show(sum(p), aspect_ratio=1)

0.000000000000000 0.100000000000000 0.200000000000000 0.300000000000000 0.400000000000000 0.500000000000000 0.600000000000000 0.700000000000000 0.800000000000000 0.900000000000000 1.00000000000000

CPU time: 0.50 s, Wall time: 0.50 s

# GSL-based solver.  This uses Runge-Kutta by default

# The ff(*yp) notation is syntax sugar for ff(yp[0], yp[1], ..., yp[n-1]) if yp has n entries.  You can think of it as splicing the list into the arguments of ff.  In our case, ff(*yp) is the same as ff(yp[0], yp[1]).

T=ode_solver()
T.function = lambda t, yp: [ff(*yp), gg(*yp)]
T.jacobian = lambda t, yp: [[fx(*yp), fy(*yp)], [gx(*yp), gy(*yp)]]

%time
T.ode_solve(y_0=initial_value, t_span=(0, t_end), num_points = num_points)
streamline = [point for t,point in T.solution]
print streamline[:10] # First 10 entries
show(line(streamline),aspect_ratio=1)

[(0, 1), [-0.0209318045202131, 0.99978090577862733], [-0.041854436965601503, 0.999123719119134], [-0.062758729280436965, 0.99802872799304498], [-0.08363552144542323, 0.99649641221281648], [-0.10447566549151191, 0.9945274432215141], [-0.12527002950843677, 0.9921226837985937], [-0.14600950164621124, 0.98928318768184209], [-0.16668499410783577, 0.98601019910564081], [-0.1872874471314655, 0.98230515225575699]]

CPU time: 0.45 s, Wall time: 0.45 s

%time
# Now we plot multiple streamlines, using GSL:
p=[]
for i in [0..1,step=0.1]:
    T.ode_solve(y_0=(0,i), t_span=(0, t_end), num_points = num_points)
    p.append(line([point for t,point in T.solution]))
show(sum(p), aspect_ratio=1)

CPU time: 0.59 s, Wall time: 0.60 s