4826-The real world problem.sagews

ACCA

Using Sage for Mathematical Modeling

Dan Warner

WWW.SAGEMATH.ORG

 

This session will provide a hands-on introduction to the basic capabilities of Sage as a tool for Mathematical Modeling

This will cover the bulk of the Maple and Mathematica capabilities that are used in Calculus and other undergraduate mathematics courses.

The real world problem

I have two sons.  Let's call them Alex and Bob.  One day I overhear the following conversation.

Bob, the younger son, said, "Give me one of your Star Wars action figures, and then we'll each have the same number."

Alex replied, "No, you give me one of yours, and I'll have twice as many as you."

I wonder how many action figures each of the boys have, and more importantly, at \$10 each, how much have I spent on this collection. 

# What don't we know?
# We don't know the number of figures that Alex has
# or the number of figures that Bob has.
# Let a denote the number of figures that Alex has
# and let b denote the number of figures that Bob has.
var('a b')

# Bob's statement can be translated as the following equation
eq1 =  a-1 == b+1

# Alex's response can be translated as the following equation
eq2 =  a+1 == 2*(b-1)

eq1;eq2

a - 1 == b + 1 a + 1 == 2*b - 2

# Solving two linear equations in two unknowns is one 
# of many straightforward algorithms that Sage implements
solve([eq1,eq2],a,b)

[[a == 7, b == 5]]

# My expense for this Star Wars habit is
10*(a+b)

10*a + 10*b

# Aha!  I need to address the fact that a and b are variables
# that would be 7 and 5 for the solution.  Here is one way to
# do that.
s = solve([eq1,eq2],a,b)
print s
10*(s[0][0]+s[0][1])

[ [a == 7, b == 5] ] 10*a + 10*b == 120

# Another way to write this is to use the dictionary form.
s = solve([eq1,eq2],a,b,solution_dict=True)
print s
10*(s[0][a]+s[0][b])

[{b: 5, a: 7}] 120

Intersecting a circle with a straight line.

Assume that the circle has radius 2 and is centered at (1,1).

Assume that the straight line passes through the points (0,-1/4) and (1/3,0).

var('x y')
eq3 = (x-1)^2 + (y-1)^2 == 4; show(eq3)
eq4 = 4*y - 3*x == -1; show(eq4)
sol3 = solve([eq3,eq4],x,y,solution_dict=True)
show(sol3)
# Use the dictionary result and display each solution to a different precision
for s in sol3:
    print(s[x].n(digits=5),s[y].n(digits=25))

\newcommand{\Bold}[1]{\mathbf{#1}}{\left(y - 1\right)}^{2} + {\left(x - 1\right)}^{2} = 4

\newcommand{\Bold}[1]{\mathbf{#1}}-3 \, x + 4 \, y = \left(-1\right)

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left\{x : -\frac{16}{25} \, \sqrt{6} + \frac{31}{25}, y : -\frac{12}{25} \, \sqrt{2} \sqrt{3} + \frac{17}{25}\right\}, \left\{x : \frac{16}{25} \, \sqrt{6} + \frac{31}{25}, y : \frac{12}{25} \, \sqrt{2} \sqrt{3} + \frac{17}{25}\right\}\right]

(-0.32767, -0.4957550765359254871346963) (2.8077, 1.855755076535925487134696)

Use the computer to display this nonlinear example and its solutions.

x0 = sol3[0][x]
y0 = sol3[0][y]
s0 = '(%5.2f,%5.2f)' % (x0,y0)
p0 = (x0-0.47,y0)
x1 = sol3[1][x]
y1 = sol3[1][y]
s1 = '(%5.2f,%5.2f)' % (x1,y1)
p1 = (x1-0.47,y1)
g = Graphics()
g += circle((1,1),2)
g += line([(-1,-1),(3,2)])
g += points([(s[x],s[y]) for s in sol3],rgbcolor=(1,0,0))
g += text(s0,p0)
g += text(s1,p1)
show(g,aspect_ratio=1)

 

# The following example by Jason Grout uses Sage to solve a system of three non-linear

# equations with 4 unknowns. We will treat u as a parameter with a value of 1

 

A nonlinear example

The following example by Jason Grout uses Sage to solve a system of three non-linear equations with 4 unknowns.

We will treat the variable $u$ as a parameter with a value of 1.

var('x y u v')
eq8 = u+v==9; show(eq8)
eq9 = v*y+u*x==-6;  show(eq9)
eq10 = v*y^2+u*x^2==24; show(eq10)

# First, we solve the system symbolically:
show(solve([eq8,eq9,eq10,u==1],u,v,x,y))

\newcommand{\Bold}[1]{\mathbf{#1}}u + v = 9

\newcommand{\Bold}[1]{\mathbf{#1}}u x + v y = \left(-6\right)

\newcommand{\Bold}[1]{\mathbf{#1}}u x^{2} + v y^{2} = 24

\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left[u = 1, v = 8, x = -\frac{4}{3} \, \sqrt{10} - \frac{2}{3}, y = \frac{1}{6} \, \sqrt{2} \sqrt{5} - \frac{2}{3}\right], \left[u = 1, v = 8, x = \frac{4}{3} \, \sqrt{10} - \frac{2}{3}, y = -\frac{1}{6} \, \sqrt{2} \sqrt{5} - \frac{2}{3}\right]\right]

# Second, we can solve it numerically
solns = solve([eq8,eq9,eq10,u==1],u,v,x,y, solution_dict=True)
for s in solns:
    print(s[u].n(digits=10), s[v].n(digits=10), s[x].n(digits=10), s[y].n(digits=10))

(1.000000000, 8.000000000, -4.883036880, -0.1396203900) (1.000000000, 8.000000000, 3.549703547, -1.193712943)

Sage does Calculus

p = (x^5-1)
diff(p^3,x)

15*(x^5 - 1)^2*x^4

integral(p^3,x)

1/16*x^16 - 3/11*x^11 + 1/2*x^6 - x

# Slightly more complicated formula using trigonometric functions with parameters
var('a b c')
q = (x-a)*(cos(b*x)+sin(c*x))
show(q)

\newcommand{\Bold}[1]{\mathbf{#1}}-{\left(\sin\left(c x\right) + \cos\left(b x\right)\right)} {\left(a - x\right)}

# Differentiate
# Illustrates the product rule, the chain rule, and the derivatives of sin and cos.
show(diff(q,x))

\newcommand{\Bold}[1]{\mathbf{#1}}{\left(a - x\right)} {\left(b \sin\left(b x\right) - c \cos\left(c x\right)\right)} + \sin\left(c x\right) + \cos\left(b x\right)

# Integrate using integration by parts
r = integral(q,x)
show(r)
show(expand(r))

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a \sin\left(b x\right)}{b} + \frac{a \cos\left(c x\right)}{c} - \frac{c x \cos\left(c x\right) - \sin\left(c x\right)}{c^{2}} + \frac{b x \sin\left(b x\right) + \cos\left(b x\right)}{b^{2}}

\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{a \sin\left(b x\right)}{b} + \frac{a \cos\left(c x\right)}{c} + \frac{x \sin\left(b x\right)}{b} - \frac{x \cos\left(c x\right)}{c} + \frac{\cos\left(b x\right)}{b^{2}} + \frac{\sin\left(c x\right)}{c^{2}}

# Define and work with functions
f(x) = sin(x)^2*cos(x)*exp(x)
show(f)

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

plot(f, 0, 3)

show(integrate(f, x))

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

x=4
show(f(x))
plot(f,0,3)

\newcommand{\Bold}[1]{\mathbf{#1}}e^{4} \sin\left(4\right)^{2} \cos\left(4\right)

A Piece of Experimental MAthematics

Let's start with a result by Archimedes, namely that

$\pi < \frac{22}{7}$.

Archimedes used a very geometric approach involving inscribed and circumscribed polygons.  

But today we might simply ask Sage, or any other modern Computer Algebra System (CAS), to evaluate the integral

$\int_0^1 \frac{(1-x)^4 x^4}{1+x^2} \,dx$

show(integral((1-x)^4*x^4/(1+x^2),x,0,1))

\newcommand{\Bold}[1]{\mathbf{#1}}-\pi + \frac{22}{7}

Since the integrand is positive on the interval (0,1), the inequality immediately follows.  

However, this may not be completely satisfactory until we look at the antiderivative.

$\int \frac{(1-x)^4 x^4}{1+x^2} \,dx$

show(integral((1-x)^4*x^4/(1+x^2),x))

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

Of course this approach not only depends on modern Computer Algebra Systems.  It also depends on the mathematical edifice of Trigonometry and Calculus that led to the wonderful relation

$\int \frac{1}{1+x^2} \,dx = arctan(x)$

$

\sum\limits_{k = 0}^\infty  {\frac{{( - 1)^{k + 1} }}{{2k + 1}}} 

$

Illustrating Partial Derivatives

var('x y u v')

(x, y, u, v)

plot3d((x^2 - y^2) / (x^2 + y^2), (x, -1, 1), (y, -1, 1))

f(x, y) = x^3 + x^2 * y^3 - 2 * y^2

P = plot3d(f(x, y), (x, 1, 3), (y, 0, 2), opacity = 0.2)
P += parametric_plot3d([2 + u, 1, f(2, 1) + diff(f(x, y), x)(x = 2, y = 1) * u], (u, -1/2, 1/2), color = "red", thickness = 2)
P += parametric_plot3d([2, 1 + u, f(2, 1) + diff(f(x, y), y)(x = 2, y = 1) * u], (u, -1/2, 1/2), color = "green", thickness = 2)
P += point((2, 1, f(2, 1)), color = "black", size = 10)
show(P)

f(x, y) = 4 - x^2 - 2 * y^2

P = parametric_plot3d([u * cos(v), 1 / sqrt(2) * u * sin(v), 4 - u^2], (u, 0, 2), (v, 0, 2 * pi), aspect_ratio = [1, 1, 1], opacity = 0.2)
P += parametric_plot3d([1 + u, 1, f(1, 1) + diff(f(x, y), x)(x = 1, y = 1) * u], (u, -1/2, 1/2), color = "red", thickness = 2)
P += parametric_plot3d([1, 1 + u, f(1, 1) + diff(f(x, y), y)(x = 1, y = 1) * u], (u, -1/2, 1/2), color = "green", thickness = 2)
P += point((1, 1, f(1, 1)), color = "black", size = 10)
show(P)

3-Dimensional Plots

Sage can draw 3d plots:

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

# Two intersecting surfaces
var('x,y')
b = 2.2
(plot3d(sin(x^2-y^2),(x,-b,b),(y,-b,b), opacity=.9) + 
    plot3d(0, (x,-b,b), (y,-b,b), color='red'))

# Implicit 3D plot
var('x,y,z')
T = golden_ratio
p = 2 - (cos(x + T*y) + cos(x - T*y) + cos(y + T*z) + cos(y - T*z) + cos(z - T*x) + cos(z + T*x))
r = 4.78
implicit_plot3d(p, (x, -r, r), (y, -r, r), (z, -r, r), plot_points=50)

Sage can plot Yoda:

from scipy import io
x = io.loadmat(DATA + 'yodapose.mat', struct_as_record=True)
from sage.plot.plot3d.index_face_set import IndexFaceSet
V = x['V']; F3 = x['F3']-1; F4 = x['F4']-1
Y = (IndexFaceSet(F3, V, color = Color('#00aa00')) + 
    IndexFaceSet(F4, V, color = Color('#00aa00')))
Y = Y.rotateX(-1)
Y.show(aspect_ratio = [1,1,1], frame = False, figsize = 4)

var('t')

parametric_plot3d([t^3, ln(3 - t), sqrt(t)], (t, 0.01, 2.99))

parametric_plot3d([cos(t), sin(t), t], (t, 0, 6 * pi))

var('s')
parametric_plot3d([cos(t), sin(t), s], (t, 0, 2 * pi), (s, 0, 6 * pi)) \
    + parametric_plot3d([cos(t), sin(t), t], (t, 0, 6 * pi), color = 'red')

parametric_plot3d([cos(t), sin(t), s], (t, 0, 2 * pi), (s, 0, 4)) \
    + parametric_plot3d([s, t, 2 - t], (s, -2, 2), (t, -2, 2), color = 'red') \
    + parametric_plot3d([cos(t), sin(t), 2 - sin(t)], (t, 0, 2 * pi), color = 'green')

parametric_plot3d([t^2, 7 * t - 12, t^2], (t, 0, 5)) \
    + parametric_plot3d([4 * t - 3, t^2, 5 * t - 6], (t, 0, 5), color = 'red')

Dynamical Systems Modeling

var('t')
y = function('y', t)
show(y)
type(y)

A Simple Growth Model

var('alpha beta')
DE1 = diff(y,t) == alpha*y - beta
show(DE1)
u1 = desolve(DE1, y, ivar=t)
show(u1)

f(t) = expand(u1)
show(f)

f(0);f(1);f(1/alpha);f(2)

g(t) = f(c=20,alpha=1/2,beta=40)
g

plot(g,0,5)

Rumor Model

Suppose that we have a somewhat isolated college town with a population of size, $P$.   At time $t=0$ a small number of people in the town learn of a rumor.  

Let $r(t)$ denote the number of people that know the rumor at time, $t$.

The spread of the rumor depends on the number of possible interactions between those who know the rumor and those who don't.  In this isolated college town the number of people who don't know the rumor at time $t$ is given by $P-r(t)$.

The possible number of interactions that result in the rumor being spread is some fraction of $r(t)*(P-r(t))$.  Let's call that fraction $\alpha$. 

Consequently, the rate of change of $r(t)$ is given by $\alpha r(t) (P-r(t))$ which gives us the following differential equation.

var('t alpha P')
r = function('r',t)
model = diff(r,t) == alpha*r*(P-r)
u3 = desolve(model,r,ivar=t)
u3

var('K')
eq = (r-P)/r == K*exp(-alpha*P*t)
eq

solve(eq,r)

We can also solve the differential equation numerically using a fourth order Runge-Kutta method.  

r=function('r',t) 
rp = desolve_rk4(diff(r,t) == (1/100)*r*(100-r), r, ics=[0,1],step=0.05, end_points=10)

A Spring Mass Damper Model

var('t k b')
x = function('x',t)
DE2 = diff(x,t,t) == -2*x - 1/2*diff(x,t) 
u2 = desolve(DE2, x, ivar=t)
u2

f(t) = u2(k1=1/2,k2=1)
f

plot(f,0,10)