3048-Sage Beginner's Guide, Chapter 1.sagews

Sage Beginner's Guide, Chapter 1

This is a cut/paste from Craig Finch's book Sage Beginner's Guide, also available on Amazon.  Amazing fun to study Sage using Sage, so to speak!

This worksheet is Chapter 1 which is an overview of many of the Sage capabilities.

x=var('x')
f(x) = (x^2 -1) / (x^4 + 1)
show(f)
show(derivative(f,x))

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{x^{2} - 1}{x^{4} + 1}

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ -\frac{4 \, {\left(x^{2} - 1\right)} x^{3}}{{\left(x^{4} + 1\right)}^{2}} + \frac{2 \, x}{x^{4} + 1}

plot(f,(-1,1))

plot(derivative(f,x),(-1,1))

f(x) = e^x * cos(x); show(f)
f_int(x) = integrate(f,x); show(f_int)

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ e^{x} \cos\left(x\right)

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{1}{2} \, {\left(\sin\left(x\right) + \cos\left(x\right)\right)} e^{x}

plot(f,(x,-2,8))

f(x) = sqrt(1 - x^2)
f_int = integrate(f, (x,0,1))
f_int.show()

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{1}{4} \, \pi

f(x) = (3 * x^4 + 4 * x^3 + 16 * x^2 + 20 * x + 9) / ((x + 2) * (x^2 + 3)^2)
g(x) = f.partial_fraction(x)
f.show();g.show()

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{3 \, x^{4} + 4 \, x^{3} + 16 \, x^{2} + 20 \, x + 9}{{\left(x + 2\right)} {\left(x^{2} + 3\right)}^{2}}

\newcommand{\Bold}[1]{\mathbf{#1}}x \ {\mapsto}\ \frac{2 \, x}{x^{2} + 3} + \frac{1}{x + 2} + \frac{4 \, x}{{\left(x^{2} + 3\right)}^{2}}

A = Matrix(QQ, [[0, -1, -1, 1], [1, 1, 1, 1], [2, 4, 1, -2], [3, 1, -2, 2]])
B = vector([0, 6, -1, 3])
show(A);show(B.column())
show(A.solve_right(B).column())

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{rrrr} 0 & -1 & -1 & 1 \\ 1 & 1 & 1 & 1 \\ 2 & 4 & 1 & -2 \\ 3 & 1 & -2 & 2 \end{array}$$\right)

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} 0 \\ 6 \\ -1 \\ 3 \end{array}$$\right)

\newcommand{\Bold}[1]{\mathbf{#1}}\left($$\begin{array}{r} 2 \\ -1 \\ 3 \\ 2 \end{array}$$\right)

f(x) = x^2
g(x) = x^3 - 2 * x^2 + 2
p1 = plot(f, (x, -1, 3), color='blue', axes_labels=['x', 'y'])
p2 = plot(g, (x, -1, 3), color='red')
show(p1+p2) # or just p1+p2

solutions=solve(f == g, x, solution_dict=True)
for s in solutions:
    show(s)

\newcommand{\Bold}[1]{\mathbf{#1}}\left\{x:\: -\sqrt{3} + 1\right\}

\newcommand{\Bold}[1]{\mathbf{#1}}\left\{x:\: \sqrt{3} + 1\right\}

\newcommand{\Bold}[1]{\mathbf{#1}}\left\{x:\: 1\right\}

labels = []; lines = []; markers = []
for s in solutions:
    x_value = s[x].n(digits=3)
    y_value = f(x_value).n(digits=3)
    labels.append(text('y=' + str(y_value), (x_value+0.5, y_value+0.5), color='black'))
    lines.append(line([(x_value, 0), (x_value, y_value)],color='black', linestyle='--'))
    markers.append(point((x_value,y_value), color='black', size=30))

show(p1+p2+sum(labels)+sum(lines)+sum(markers))

var('u,v')
r = 2.0
f_x = (r + cos(u / 2) * sin(v) - sin(u / 2) * sin(2 * v)) * cos(u)
f_y = (r + cos(u / 2) * sin(v) - sin(u / 2) * sin(2 * v)) * sin(u)
f_z = sin(u / 2) * sin(v) + cos(u / 2) * sin(2 * v)
parametric_plot3d([f_x, f_y, f_z], (u, 0, 2 * pi), (v, 0, 2 * pi), color="red")

import numpy
var('OD, slope, intercept')
def standard_curve(OD, slope, intercept): 
    """Apply a linear standard curve to optical density data""" 
    return OD * slope + intercept
# Enter data to define standard curve 
CFU = numpy.array([2.60E+08, 3.14E+08, 3.70E+08, 4.62E+08, 8.56E+08, 1.39E+09, 1.84E+09]) 
optical_density = numpy.array([0.083, 0.125, 0.213, 0.234, 0.604, 1.092, 1.141]) 
OD_vs_CFU = zip(optical_density, CFU)
# Fit linear standard 
std_params = find_fit(OD_vs_CFU, standard_curve,
    parameters=[slope, intercept], 
    variables=[OD], initial_guess=[1e9, 3e8], 
    solution_dict = True)
for param, value in std_params.iteritems():
    print(str(param) + ' = %e' % value)
# Plot 
data_plot = scatter_plot(OD_vs_CFU, markersize=20,
    facecolor='red', axes_labels=['OD at 600nm', 'CFU/ml'])
fit_plot = plot(standard_curve(OD, std_params[slope], std_params[intercept]), (OD, 0, 1.2))
show(data_plot+fit_plot)

intercept = 1.237283e+08 slope = 1.324714e+09

sample_times = numpy.array([0, 20, 40, 60, 80, 100, 120, 140, 160, 180, 200, 220, 240, 280, 360, 380, 400, 420, 440, 460, 500, 520, 540, 560, 580, 600, 620, 640, 660, 680, 700, 720, 760, 1240, 1440, 1460, 1500, 1560])
OD_readings = numpy.array([0.083, 0.087, 0.116, 0.119, 0.122, 0.123, 0.125, 0.131, 0.138, 0.142, 0.158, 0.177, 0.213, 0.234, 0.424, 0.604, 0.674, 0.726, 0.758, 0.828, 0.919, 0.996, 1.024, 1.066, 1.092, 1.107, 1.113, 1.116, 1.12, 1.129, 1.132, 1.135, 1.141, 1.109, 1.004, 0.984, 0.972, 0.952])
concentrations = standard_curve(OD_readings, std_params[slope], std_params[intercept])
exp_data = zip(sample_times, concentrations)
data_plot = scatter_plot(exp_data, markersize=20, facecolor='red', axes_labels=['time (sec)', 'CFU/ml'])
show(data_plot)

var('t, max_rate, lag_time, y_max, y0')
def gompertz(t, max_rate, lag_time, y_max, y0): 
    """Define a growth model based upon the Gompertz growth curve""" 
    return y0 + (y_max - y0) * numpy.exp(-numpy.exp(1.0 + 
    max_rate * numpy.exp(1) * (lag_time - t) / (y_max - y0)))
# Estimated parameter values for initial guess 
max_rate_est = (1.4e9 - 5e8)/200.0 
lag_time_est = 100 
y_max_est = 1.7e9
y0_est = 2e8
gompertz_params = find_fit(exp_data, gompertz, 
    parameters=[max_rate, lag_time, y_max, y0], 
    variables=[t], 
    initial_guess=[max_rate_est, lag_time_est, y_max_est, y0_est], 
    solution_dict = True)
for param,value in gompertz_params.iteritems(): 
    print(str(param) + ' = %e' % value)

y0 = 3.051120e+08 y_max = 1.563822e+09 lag_time = 2.950639e+02 max_rate = 6.606311e+06

gompertz_model_plot = plot(gompertz(t, gompertz_params[max_rate], 
    gompertz_params[lag_time], gompertz_params[y_max], 
    gompertz_params[y0]), (t, 0, sample_times.max()))
show(gompertz_model_plot + data_plot)