# For an introduction to Sage work through the worksheet, CE_Math_Intro_01 # # 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, e.g. 'plot?' tells you about plot 

# 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 CAS. 

# 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 (pardon the pun) 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) 

# 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) 

# You can substitute values for all or just some parameters show(q0(a=-1,c=1)) show(q0.subs(b=5)) # These expressions can be plotted 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 does know 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) 

# Surface plot example var('x,y') plot3d(cos(x^2 + y^2), (x, -2, 2), (y, -2, 2)) 

# 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]) 

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

var('s') print x(s) expand(x(s));x(0);x(0)==4;solve(x(0)==4,c) 

# Working with Ordinary Differential Equations var('u t') ODE(t, u) = - 1.5 * (u - 60 - 15 * sin(2 * pi * t)) VF = plot_slope_field(ODE(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: [ODE(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) 

# Plot a vector field for the phase plane portrait # for a simple predator prey model var('r w') ODEr(r, w) = r - w * r ODEw(r, w) = 0.25 * w * r - w VF = plot_vector_field([ODEr(r, w), ODEw(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: [ODEr(y[0], y[1]), ODEw(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')