Computational and Experimental Mathematics

Calculus with SAGE

Neil Calkin, Dan Warner, and Holly Hirst

November 14, 2010

# For an introduction to Sage work through the worksheet, SC10CE1a # # The Sage notebook uses blocks of text, which are essentially Python statements that get executed. # Lines that start with the pound sign are comments and do not get echoed when the block is evaluated. # # A block is evaluated by clicking on the evaluate link or by typing shift-enter. # # Pressing tab will give auto-complete options, or function descriptions after a parens ( # A '?' after a command gives detailed information, for example, 'plot?' tells you about the plot command.

###

Sage can do Calculus

p0 = (x^2 - 1/2)^3 show(p0) p1 = diff(p0,x) show(p1) p2 = integral(p0,x) show(p2)

Hard nosed Calculus teachers would not accept this last answer,
but leaving off the 'c' is common in Computer Algebra Systems.

###

These expressions can be plotted.

P0 = plot( p0, -0.5, 0.5, rgbcolor=(1,0,0) ) # red P1 = plot( p1, -0.5, 0.5, rgbcolor=(0,1,0) ) # green P2 = plot( p2, -0.5, 0.5, rgbcolor=(0,0,1) ) # blue show(P0+P1+P2)

###

Sage handles rational functions with ease

r0 = x / (x^2 - 1) show(r0) r1 = diff(r0,x) show(r1) r2 = diff(r0,x,2) show(r2)

These expressions can also be plotted.

P0 = plot( r0, -0.5, 0.5, rgbcolor=(1,0,0) ) P1 = plot( r1, -0.5, 0.5, rgbcolor=(0,1,0) ) P2 = plot( r2, -0.5, 0.5, rgbcolor=(0,0,1) ) show(P0+P1+P2)

Sage can make substitutions, generate partial fraction expansions, and integrate.

var('z') r3 = r0(x=z) show(r3) pf = r3.partial_fraction(z) show(pf) r4 = integral(pf,z) show(r4)

Most Calculus teachers would not accept this last answer, even though it is
correct if log is the natural log and the antiderivative is defined for complex z.
Nonetheless, the lack of context can lead to real (pun intended) problems.

###

On the other hand, Sage can nicely handle this important integral

show(4*integral(1/(1+x^2),x)) show(4*integral(1/(1+x^2),x,0,1)) print pi.n(digits=140)

As an aside, here is a fun web site that searches the first 200 million digits of $\pi$ for requested strings of digits, such as your birthday.

http://www.angio.net/pi/bigpi.cgi

It provides a great take-home activity for elementary talk about the history of computing $\pi$.

Here's a slightly more complicated formula with trigonometric functions with parameters.

var('a b c') q0 = (x-a)*(cos(b*x)+sin(c*x)) show(q0) print type(q0)

Differentiate using the product rule, the chain rule, and the derivatives of sin and cos.
Note the variety of ways to tell Sage to differentiate.

show(diff(q0,x)) show(q0.diff(x)) # diff is both a function and a method show(derivative(q0,x)) show(q0.derivative(x)) # derivative is both a function and a method show(q0.differentiate(x)) # differentiate is only a method # show(differentiate(q0,x)) q1 = diff(q0,x)

Integrate using integration by parts.

q2 = integral(q0,x) show(q2) show(expand(q2))

The expanded form is much closer to the typical answer done by hand.

You can substitute values for all or just some parameters.

show(q0(a=-1,c=1)) show(q0.subs(b=5))

Of course, these expressions can be plotted when all the parameters are specified.

qq0 = q0(a=-1,b=2,c=1) qq1 = q1(a=-1,b=2,c=1) qq2 = q2(a=-1,b=2,c=1) P0 = plot( qq0, -1.0, 1.0, rgbcolor=(1,0,0) ) P1 = plot( qq1, -1.0, 1.0, rgbcolor=(0,1,0) ) P2 = plot( qq2, -1.0, 1.0, rgbcolor=(0,0,1) ) show(P0+P1+P2)

Sage knows about the error function.

f0 = exp(-x^2) show(f0) f1 = integral(f0,x) show(f1) P0 = plot(f0, -2.0, 2.0, rgbcolor=(1,0,0)) P1 = plot(f1, -2.0, 2.0, rgbcolor=(0,0,1)) show(P0+P1)

###

Sage can generate surface plots.

var('x,y') plot3d(cos(x^2 + y^2), (x, -2, 2), (y, -2, 2))

Sage can generate parametric plots.

# Set up the parameters var('u v') # The x, y, & z coordinates as functions of the parameters fx = cos(u) * sin(v) fy = 2 * sin(u) * sin(v) fz = 3 * cos(v) # Plot an ellipsoid # 'aspect_ratio = [1, 1, 1]' makes the scaling of the axes equal # without it, the plot looks like a sphere because the axes are each scaled differently parametric_plot3d([fx, fy, fz], (u, 0, 2 * pi), (v, 0, pi), aspect_ratio = [1, 1, 1])

###

Sage can symbolically solve ODEs

Consider the Initial Value Problem

$ \frac{dx(t)}{dt}= 1 - x(t)$ with $x(0) = 4$

t = var('t') # define a variable t x = function('x',t) # define x to be a function of that variable print x DE1 = diff(x, t) == 1 - x # define the ordinary differential equation x(t) = desolve(DE1, x) # Solve the ODE symbolically using Maxima print 'Solution: ',x

# Sage can also provide feedback on the method used. u = function('u',t) u(t) = desolve(diff(u,t)==1-u, u, show_method=true) print 'Solution: ',u

var('c') print x(t) print expand(x(t)) # This is the form we would customarily expect. print x(0) print x(0)==4 # Write an equation for an initial condition solve(x(0)==4,c) # Solve the initial condition

x = function('x',t) x(t) = desolve(diff(x,t)==1-x, [x,t], ics=[0,4]) # Solve the IVP symbolically using Maxima print 'Solution: ',x

Solve the Rumor Model - a first order nonlinear ODE.

var('t P r') # define a variable t x = function('x',t) # define x to be a function of that variable RM = diff(x,t) == r*x*(P-x) x(t) = desolve(RM, [x, t]) # Solve the ODE symbolically using Maxima print 'Solution: ',x

Solve a second order linear ODE symbolically.

t = var('t') # define a variable t x = function('x',t) # define x to be a function of that variable DE = diff(x, t, 2) - x == t # define the ordinary differential equation x(t) = desolve(DE, x) # Solve the ODE symbolically using Maxima print 'Solution: ',x

x = function('x',t) x(t) = desolve(DE, [x,t], ics=[10,2,1]) # Solve the IVP symbolically using Maxima print 'Solution: ',x

###

Sage can solve ODE's numerically

It provides two different packages - the desolve library and the ode_solver class.

# desolve with only the right hand side specified var('x t') rhs(t,x) = 1-x desolve_rk4(rhs, x, ivar=t, ics=[0,4], end_points=1, step=0.25)

# desolve with the ODE specified x = function('x',t) DE1 = diff(x, t) == 1 - x desolve_rk4(DE1, x, ics=[0,4], end_points=1, step=0.25)

# Using the ode_solve class T = ode_solver() var('t x') rhs(t,x) = 1 - x T.function = lambda t, x: [rhs(t, x[0])] T.ode_solve(y_0 = [4], t_span = [0, 1], num_points = 4) T.solution

###

Here are some common manipulations for working with ODE's numerically.

var('u t') rhs(t, u) = - 3/2 * (u - 60 - 15 * sin(2 * pi * t)) print type(rhs) show(rhs)

Plot the slope field.

VF = plot_slope_field(rhs(t, u), [t, 0, 3], [u, 0, 90]) show(VF)

Compute several solutions and display them on the slope field.

T = ode_solver() T.function = lambda t, u: [rhs(t, u[0])] # y_0 = 10 T.ode_solve(y_0 = [10], t_span = [0, 3], num_points = 5) T.solution

T = ode_solver() T.function = lambda t, u: [rhs(t, u[0])] solutions = [] for y_0 in range(0, 100, 10): T.ode_solve(y_0 = [y_0], t_span = [0, 3], num_points = 100) solutions.append(line([[p[0], p[1][0]] for p in T.solution])) show(sum(solutions) + VF)

###

First order systems of ODEs

Plot a vector field of the phase plane portrait
for a simple predator prey model

var('r w') ppr(r, w) = r - w * r ppw(r, w) = 0.25 * w * r - w VF = plot_vector_field([ppr(r, w), ppw(r, w)], [r, 0, 10], [w, 0, 3]) show(VF)

Compute several solutions and display them on the vector field

T = ode_solver() T.function = lambda t, y: [ppr(y[0], y[1]), ppw(y[0], y[1])] solutions = [] for y_0 in range(4, 11): T.ode_solve(y_0 = [y_0, 1], t_span = [0, 10], num_points = 100) solutions.append(line([p[1] for p in T.solution])) equilib = [4, 1] show(sum(solutions) + VF)

Calculate one solution set and display it on a standard time plot.

T.ode_solve(y_0 = [7, 1], t_span = [0, 20], num_points = 1000) line([[p[0], p[1][0]] for p in T.solution]) + line([[p[0], p[1][1]] for p in T.solution], rgbcolor = 'red')

SC10CE2a

4583 days ago by SC2010