CoCalc
Sharedscratch.sagewsOpen in CoCalc
Author: Paul Bryan


+_
from sage.plot.plot3d.implicit_surface import ImplicitSurface
var('x,y,z')
t = (1-sin(2*x*y+3*z)**2).function(x,y,z)
cm = colormaps.autumn
G = ImplicitSurface(x^2 + y^2 + z^2, (x,-2, 2), (y,-2, 2), (z,-2, 2), contour=4, color=(t,cm))
G.show()
G.show(viewer='tachyon')

(x, y, z)
3D rendering not yet implemented
f(x) = x*exp(x)/(exp(x) - 1)
plot(f, (x, 0, 1))

# Solve Riccati equation
t = var('t')
A = var('A')
p = var('p')
q = function('q')(t)
L, K, B = var('L, K, B')

# Flat space
ode = diff(q, t) == (p+1)/p * q^2
ode_soln = desolve(ode, q, ivar=t)

show("Flat solution")
show(ode)
#show(ode_soln)
soln = solve(ode_soln, q)
show(soln[0])

latex(soln[0])

show(soln[0].substitute(p = 1))
latex(soln[0].substitute(p = 1))

Flat solution
$\displaystyle \frac{\partial}{\partial t}q\left(t\right) = \frac{{\left(p + 1\right)} q\left(t\right)^{2}}{p}$
$\displaystyle q\left(t\right) = -\frac{p}{C {\left(p + 1\right)} + {\left(p + 1\right)} t}$
q\left(t\right) = -\frac{p}{C {\left(p + 1\right)} + {\left(p + 1\right)} t}
$\displaystyle q\left(t\right) = -\frac{1}{2 \, {\left(C + t\right)}}$
q\left(t\right) = -\frac{1}{2 \, {\left(C + t\right)}}

# Sphere
#ode = diff(q, t) == (p+1)/p * q^2 -(4*A/p) * q + (2/p) * A^2
ode = diff(q, t) == L * q^2 + K * q + B

assume(4 * B * L - K^2 < 0)
ode_soln = desolve(ode, q, ivar=t)

show("Sphere solution")
show(ode)
show(ode.substitute(L = (p+1)/p, K = - (4/p) * A, B = (2/p) * A^2))
#show(ode_soln)
soln = solve(ode_soln, q)
show(soln[0].substitute(C = 0))
show(soln[0].substitute(L = (p+1)/p, K = - (4/p) * A, B = (2/p) * A^2))

latex(soln[0])

Sphere solution
$\displaystyle \frac{\partial}{\partial t}q\left(t\right) = L q\left(t\right)^{2} + K q\left(t\right) + B$
$\displaystyle \frac{\partial}{\partial t}q\left(t\right) = \frac{{\left(p + 1\right)} q\left(t\right)^{2}}{p} + \frac{2 \, A^{2}}{p} - \frac{4 \, A q\left(t\right)}{p}$
$\displaystyle q\left(t\right) = -\frac{{\left(K + \sqrt{K^{2} - 4 \, B L}\right)} e^{\left(\sqrt{K^{2} - 4 \, B L} C + \sqrt{K^{2} - 4 \, B L} t\right)} - K + \sqrt{K^{2} - 4 \, B L}}{2 \, {\left(L e^{\left(\sqrt{K^{2} - 4 \, B L} C + \sqrt{K^{2} - 4 \, B L} t\right)} - L\right)}}$
$\displaystyle q\left(t\right) = -\frac{{\left(\sqrt{-\frac{8 \, A^{2} {\left(p + 1\right)}}{p^{2}} + \frac{16 \, A^{2}}{p^{2}}} - \frac{4 \, A}{p}\right)} e^{\left(C \sqrt{-\frac{8 \, A^{2} {\left(p + 1\right)}}{p^{2}} + \frac{16 \, A^{2}}{p^{2}}} + \sqrt{-\frac{8 \, A^{2} {\left(p + 1\right)}}{p^{2}} + \frac{16 \, A^{2}}{p^{2}}} t\right)} + \sqrt{-\frac{8 \, A^{2} {\left(p + 1\right)}}{p^{2}} + \frac{16 \, A^{2}}{p^{2}}} + \frac{4 \, A}{p}}{2 \, {\left(\frac{{\left(p + 1\right)} e^{\left(C \sqrt{-\frac{8 \, A^{2} {\left(p + 1\right)}}{p^{2}} + \frac{16 \, A^{2}}{p^{2}}} + \sqrt{-\frac{8 \, A^{2} {\left(p + 1\right)}}{p^{2}} + \frac{16 \, A^{2}}{p^{2}}} t\right)}}{p} - \frac{p + 1}{p}\right)}}$
q\left(t\right) = -\frac{{\left(K + \sqrt{K^{2} - 4 \, B L}\right)} e^{\left(\sqrt{K^{2} - 4 \, B L} C + \sqrt{K^{2} - 4 \, B L} t\right)} - K + \sqrt{K^{2} - 4 \, B L}}{2 \, {\left(L e^{\left(\sqrt{K^{2} - 4 \, B L} C + \sqrt{K^{2} - 4 \, B L} t\right)} - L\right)}}
# Test analytic solution
a = -(2/p)*A
b = sqrt((A)/p - (p+1)*A^2/(2*p))

a = -K/2
b = (1/2) * sqrt(K^2 - 4 * B *L)

q_anal(t) = (1/L) * (a + b - (a-b)*exp(2*b*t))/(1 - exp(2*b*t))

q1(t) = soln[0].rhs()

show("Analytic Solution")
show(q_anal)

show(ode)
#ode.substitute_function(q, q0)
(q_anal.diff(t) - L * q_anal^2 - K * q_anal - B).full_simplify()

show("Taylor Series")
q_taylor = (t * q_anal).substitute(L = (p+1)/p, K = - (4/p) * A, B = (2/p) * A^2)
show(Yq_taylor.taylor(t, 0, 3))

Analytic Solution
$\displaystyle t \ {\mapsto}\ -\frac{{\left(K + \sqrt{K^{2} - 4 \, B L}\right)} e^{\left(\sqrt{K^{2} - 4 \, B L} t\right)} - K + \sqrt{K^{2} - 4 \, B L}}{2 \, L {\left(e^{\left(\sqrt{K^{2} - 4 \, B L} t\right)} - 1\right)}}$
$\displaystyle \frac{\partial}{\partial t}q\left(t\right) = L q\left(t\right)^{2} + K q\left(t\right) + B$
0
Taylor Series
$\displaystyle t \ {\mapsto}\ \frac{2 \, A^{2} {\left(p - 1\right)} t^{2}}{3 \, {\left(p^{2} + p\right)}} + \frac{2 \, A t}{p + 1} - \frac{p}{p + 1}$
# Plot the solution
p0 = 1/2
A0 = 3
L0 = (p0+1/p0)
K0 = - (4/p0) * A0
B0 = (2/p0) * A0^2

show(L0)
show(K0)
show(B0)
show(K0^2 - 4 * B0 * L0)

q0 = q_anal.substitute(L = L0, K = K0, B = B0)

show(q0)

plot(t*q0, (t, 0, 1))


$\displaystyle \frac{5}{2}$
$\displaystyle -24$
$\displaystyle 36$
$\displaystyle 216$
$\displaystyle t \ {\mapsto}\ -\frac{6 \, {\left({\left(\sqrt{6} - 4\right)} e^{\left(6 \, \sqrt{6} t\right)} + \sqrt{6} + 4\right)}}{5 \, {\left(e^{\left(6 \, \sqrt{6} t\right)} - 1\right)}}$

# Sphere throwing away A^2 > 0 term
#ode = diff(q, t) == (p+1)/p * q^2 -(4*A/p) * q
ode = diff(q, t) == C * q^2 + K * q
assume(p-1.0>0)
ode_soln = desolve(ode, q, ivar=t)

show("Sphere solution without $A^2$")
show(ode)
show(ode_soln)
soln = solve(ode_soln, q)[0]
show(soln)
show(simplify(soln.rhs()))

Sphere solution without $A^2$
$\displaystyle \frac{\partial}{\partial t}q\left(t\right) = C q\left(t\right)^{2} + K q\left(t\right)$
$\displaystyle -\frac{\log\left(C q\left(t\right) + K\right) - \log\left(q\left(t\right)\right)}{K} = C + t$
$\displaystyle \log\left(q\left(t\right)\right) = K C + K t + \log\left(C q\left(t\right) + K\right)$
$\displaystyle K C + K t + \log\left(C q\left(t\right) + K\right)$

# Test analytic solution wihtout A^2
q0(t) = (K/C) * exp(K*t)/(1 - exp(K*t))
show(q0)
show(ode)
#ode.substitute_function(q, q0)

q0.diff(t) - C * q0^2 - K * q0

$\displaystyle t \ {\mapsto}\ -\frac{K e^{\left(K t\right)}}{C {\left(e^{\left(K t\right)} - 1\right)}}$
$\displaystyle \frac{\partial}{\partial t}q\left(t\right) = C q\left(t\right)^{2} + K q\left(t\right)$
t |--> 0


p = plot([])
p += line([[-1,0], [1, 0]])
p += point((-1,0), size=50)
p += point((1,0), size=50)

p.show(axes=False)
p.save("straight.png", axes=False)

p = plot([])
p += arc((0,0), 1, sector=(0, pi))
p += line([[-1,0], [1, 0]])
p += point((-1,0), size=50)
p += point((1,0), size=50)

p.show(axes=False)
p.save("large_curvature.png", axes=False)

p = plot([])
r = 3
p += arc((0,-sqrt(r^2-1)), r, sector=(pi/2-arcsin(1/r), pi/2+arcsin(1/r)))
p += line([[-1,0], [1, 0]])
p += point((-1,0), size=50)
p += point((1,0), size=50)

p.show(axes=False, ymin=-0.1)
p.save("small_curvature.png", axes=False, ymin=-0.1)


p = plot([])
p += bezier_path([[(0,0), (1,0), (0, -1), (1,0)]], color="blue")

p += line([[0,0], [1, 0]])
p += point((0,0), size=50)
p += point((1,0), size=50)

p.show(axes=False)

integrate(1/sin(x), x)

-1/2*log(cos(x) + 1) + 1/2*log(cos(x) - 1)


u(x) = -(1/2) * x + x^2
up(x) = u.diff(x)

p = u - up

show(p)
show((5/2)^2 - 2)
x0 = find_root(p, 0, 10)

q = plot(up - u, (x, 0, 10))
q += point((x0, 0))
q.show()


$\displaystyle x \ {\mapsto}\ x^{2} - \frac{5}{2} \, x + \frac{1}{2}$
$\displaystyle \frac{17}{4}$





(x, t, u) = var('x,t, u')
k(x, t) = sqrt((1 - exp(2*t))^(-1) - (sin(x))^2)

assume(0 < u < pi/2)
assume(t < 0)
assume(exp(t)-1>0)
y(u, t) = integrate(cos(x)/k(x, t), (x, 0, u))
z(u, t) = integrate(sin(x)/k(x, t), (x, 0, u))

show(k)
show(y)
show(z)


$\displaystyle \left( x, t \right) \ {\mapsto} \ \sqrt{-\sin\left(x\right)^{2} - \frac{1}{e^{\left(2 \, t\right)} - 1}}$
$\displaystyle \left( u, t \right) \ {\mapsto} \ -i \, {\rm arcsinh}\left(\frac{\sin\left(u\right)}{\frac{1}{\sqrt{e^{\left(2 \, t\right)} - 1}}}\right)$
$\displaystyle \left( u, t \right) \ {\mapsto} \ -\log\left(2 \, \sqrt{\frac{{\left(e^{\left(2 \, t\right)} - 1\right)} \cos\left(u\right)^{2} - e^{\left(2 \, t\right)}}{e^{\left(2 \, t\right)} - 1}} + 2 \, \cos\left(u\right)\right) + \log\left(2 \, \sqrt{-\frac{1}{e^{\left(2 \, t\right)} - 1}} + 2\right)$
# Plot integer lattice
M = 4
w = 1/4

lit_bg = "aqua"
unlit_bg = "lightgray"
lit_fg = "black"
unlit_fg = "gray"

p = plot([])

lattice = [(m, n) for m in range(-M, M+1) for n in range(-M, M+1)]

def lightup(m, n):
if x == 0 and y == 0: return True
if y < 0: return False
if gcd(x, y) == 1: return True

for l in lattice:
x = l[0]
y = l[1]
p += polygon([(x-w,y-w), (x-w,y+w), (x+w, y+w), (x+w, y-w)], color = lit_bg if lightup(x, y) else unlit_bg)
p += text("%d/%d" % (x, y), (x, y), color = lit_fg if (gcd(x, y) == 1 or (x == 0 and y == 0)) else unlit_fg)

p.show(aspect_ratio=1, axes=False)

import itertools

arr = [[1,2,3,4],
[12,13,14,5],
[11,16,15,6],
[10,9,8,7]]

arr = lattice

def transpose_and_yield_top(arr):
while arr:
yield arr[0]
arr = list(reversed(zip(*arr[1:])))

for a in arr:
print(a)

l = list(itertools.chain(*transpose_and_yield_top(arr)))
#l.reverse()
print l


(-4, -4) (-4, -3) (-4, -2) (-4, -1) (-4, 0) (-4, 1) (-4, 2) (-4, 3) (-4, 4) (-3, -4) (-3, -3) (-3, -2) (-3, -1) (-3, 0) (-3, 1) (-3, 2) (-3, 3) (-3, 4) (-2, -4) (-2, -3) (-2, -2) (-2, -1) (-2, 0) (-2, 1) (-2, 2) (-2, 3) (-2, 4) (-1, -4) (-1, -3) (-1, -2) (-1, -1) (-1, 0) (-1, 1) (-1, 2) (-1, 3) (-1, 4) (0, -4) (0, -3) (0, -2) (0, -1) (0, 0) (0, 1) (0, 2) (0, 3) (0, 4) (1, -4) (1, -3) (1, -2) (1, -1) (1, 0) (1, 1) (1, 2) (1, 3) (1, 4) (2, -4) (2, -3) (2, -2) (2, -1) (2, 0) (2, 1) (2, 2) (2, 3) (2, 4) (3, -4) (3, -3) (3, -2) (3, -1) (3, 0) (3, 1) (3, 2) (3, 3) (3, 4) (4, -4) (4, -3) (4, -2) (4, -1) (4, 0) (4, 1) (4, 2) (4, 3) (4, 4) [-4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, -4, -3, -2, -1, 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -2, -2, -2, -2, -2, -2, -2, -2, -2, -3, -3, -3, -3, -3, -3, -3, -3, -3, -4, -4, -4, -4, -4, -4, -4, -4]


from fipy import *

nx = 50
dx = 1

mesh = Grid1D(nx = nx, dx = dx)


Error in lines 4-4 Traceback (most recent call last): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 995, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/fipy/meshes/factoryMeshes.py", line 187, in Grid1D if numerix.getShape(dx) == (): File "/projects/sage/sage-7.5/local/lib/python2.7/site-packages/fipy/tools/numerix.py", line 242, in getShape raise AttributeError, "No attribute 'shape'" AttributeError: No attribute 'shape'
N = 200
W = 1
H = 1
seglength = 5

curve = []
for i in range(N):
ci = cos((2*pi)*(i/N))
si = sin((2*pi)*(i/N))

u = W/2 + W * (0.05 * ci)
x = u + 0.2 * W * ci^101
y = H * (0.15 + 0.05 * si + 0.05 * sin(u/5) + 0.7 * si^150)

curve += [(x, y)]

z = var('z')
f(z) = W/2 + W * 0.05 * cos(z)
analytic_curve(z) = (f(z) + 0.2 * W * cos(z)^101, H * (0.15 + 0.05 * sin(z) + 0.05 * sin(f(z)/5) + 0.7 * sin(z)^150))

for i in range(len(curve)):
a = curve[i]
b = curve[i+1 if i < len(curve) - 1 else 0]

dx = b[0] - a[0]
dy = b[1] - a[1]

dr2 = dx^2 + dy^2

if dr2 > 4 * seglength^2:
dr = sqrt(dr2)
else:
if len(curve) > 4 and 2 * dr2 < seglength^2:
dr = sqrt(dr2)

#var x = canvas.width/2 + canvas.width*(0.05 * Math.cos(2*Math.PI*i/N));
#        curve.push({
#            x: x + 0.2*canvas.width*Math.pow(Math.cos(2*Math.PI*i/N),101),
#            y: canvas.height * (0.15 + 0.05 * Math.sin(2*Math.PI*i/N) + 0.05*Math.sin(x/5) + 0.7 * Math.pow(Math.sin(2*Math.PI*i/N), 150))
#        });

#for (var i = 0; i < cu.length; i++) {
#        var a = cu[i];
#        var bi = (i < cu.length - 1 ? i+1 : 0), b = cu[bi];
#
#        var dx = b.x - a.x;
#        var dy = b.y - a.y;
#
#        var dr2 = dx*dx + dy*dy;
#        if (dr2 > 4*seglength*seglength) {
#            var dr = Math.pow(dr2, 1/2);
#            cu.splice(1+i,0,{
#                x: a.x + seglength * dx/dr,
#                y: a.y + seglength * dy/dr
#            });
#        }
#        else if (cu.length > 4 && dr2 * 2 < seglength * seglength) {
#            cu.splice(i--,1);
#        }
#    }

p = plot([])

for i in range(len(curve)):
j = i+1 if i < len(curve) - 1 else 0
p += line([curve[j-1], curve[j]])

p.show()

q = parametric_plot(analytic_curve(z), (z, 0, 2*pi))
q.show(aspect_ratio=1)


c = var('c')
G(x) = ((sin(x))^2 + c)^(-1/2)
G(x) = (sin(x))^(-1)
H(x) = integrate(G, x)

show(G)
show(H)

$\displaystyle x \ {\mapsto}\ \frac{1}{\sin\left(x\right)}$
$\displaystyle x \ {\mapsto}\ -\frac{1}{2} \, \log\left(\cos\left(x\right) + 1\right) + \frac{1}{2} \, \log\left(\cos\left(x\right) - 1\right)$
import wave
from sage.media.wav import Wave
w = Wave('Alesis-Fusion-Voice-Oohs-C4.wav')

w.listen()
p = w.plot()
p.save('wav.png')
p = w.plot_fft()
p.save('wav_ft.png')

Click to listen to Alesis-Fusion-Voice-Oohs-C4.wav
# Plot a cylinder
from sage.plot.plot3d.shapes import Cylinder

(x,y,z) = var('x,y,z')

p = plot([])

t = (z).function(x,y,z)
cmap = colormaps.coolwarm
#p += ImplicitSurface(x^2 + y^2, (x,-1, 1), (y,-1, 1), (z, 0, 1), contour=1/4, color=(t,cm))
p += arrow3d((1, 0, 0), (1, 0, 1/2), color="black")

p += text3d("q", (1.2,0,0), fontweight="bold")

#p.show(frame=False)
#p.save("fourier_law.png", frame=False)

p = plot([])

t = (y).function(x, y, z)
p += contour_plot(y, (x, 0, 1/2), (y, 0, 1), cmap=cmap, contours=500)

p += text(r"$\vec{q} = - \nabla T$", (3/4, 1/2), fontweight="bold", fontsize="large", color="black")
p += line([(0, 1/2), (1/2, 1/2)], color="black", linestyle="--")
for i in range(5):
p += arrow((i/8, 1/2), (i/8, 3/8), color="black", width=1, arrowsize=3)

p.show(frame=False)
p.save("fourier_law.png", frame=False)


sorted(colormaps)

arrow3d.?


[u'Accent', u'Accent_r', u'Blues', u'Blues_r', u'BrBG', u'BrBG_r', u'BuGn', u'BuGn_r', u'BuPu', u'BuPu_r', u'CMRmap', u'CMRmap_r', u'Dark2', u'Dark2_r', u'GnBu', u'GnBu_r', u'Greens', u'Greens_r', u'Greys', u'Greys_r', u'OrRd', u'OrRd_r', u'Oranges', u'Oranges_r', u'PRGn', u'PRGn_r', u'Paired', u'Paired_r', u'Pastel1', u'Pastel1_r', u'Pastel2', u'Pastel2_r', u'PiYG', u'PiYG_r', u'PuBu', u'PuBuGn', u'PuBuGn_r', u'PuBu_r', u'PuOr', u'PuOr_r', u'PuRd', u'PuRd_r', u'Purples', u'Purples_r', u'RdBu', u'RdBu_r', u'RdGy', u'RdGy_r', u'RdPu', u'RdPu_r', u'RdYlBu', u'RdYlBu_r', u'RdYlGn', u'RdYlGn_r', u'Reds', u'Reds_r', u'Set1', u'Set1_r', u'Set2', u'Set2_r', u'Set3', u'Set3_r', u'Spectral', u'Spectral_r', u'Wistia', u'Wistia_r', u'YlGn', u'YlGnBu', u'YlGnBu_r', u'YlGn_r', u'YlOrBr', u'YlOrBr_r', u'YlOrRd', u'YlOrRd_r', u'afmhot', u'afmhot_r', u'autumn', u'autumn_r', u'binary', u'binary_r', u'bone', u'bone_r', u'brg', u'brg_r', u'bwr', u'bwr_r', u'cool', u'cool_r', u'coolwarm', u'coolwarm_r', u'copper', u'copper_r', u'cubehelix', u'cubehelix_r', u'flag', u'flag_r', u'gist_earth', u'gist_earth_r', u'gist_gray', u'gist_gray_r', u'gist_heat', u'gist_heat_r', u'gist_ncar', u'gist_ncar_r', u'gist_rainbow', u'gist_rainbow_r', u'gist_stern', u'gist_stern_r', u'gist_yarg', u'gist_yarg_r', u'gnuplot', u'gnuplot2', u'gnuplot2_r', u'gnuplot_r', u'gray', u'gray_r', u'hot', u'hot_r', u'hsv', u'hsv_r', u'jet', u'jet_r', u'nipy_spectral', u'nipy_spectral_r', u'ocean', u'ocean_r', u'pink', u'pink_r', u'prism', u'prism_r', u'rainbow', u'rainbow_r', u'seismic', u'seismic_r', u'spectral', u'spectral_r', u'spring', u'spring_r', u'summer', u'summer_r', u'terrain', u'terrain_r', u'winter', u'winter_r']
File:
Docstring :


theta = var('theta')

x(theta) = cos(theta) + 1.2
y(theta) = sin(theta) - 1/2

parametric_plot((x, y), (theta, 0, 2*pi))

var('s,t')
def sphere_and_plane(x):
return sphere((0,0,0),1,color='red',opacity=.5)+parametric_plot3d([t,x,s],(s,-1,1),(t,-1,1),color='green',opacity=.7)

sp = animate([sphere_and_plane(x) for x in sxrange(-1,1,.3)], frame=False)
sp[0]      # first frame
sp[-1]     # last frame

sp.show()  # optional -- ImageMagick


(s, t)
3D rendering not yet implemented
3D rendering not yet implemented

(x,y,z) = var('x,y,z')
def frame(t):
return implicit_plot3d((x^2 + y^2 + z^2), (x, -2, 2), (y, -2, 2), (z, -2, 2), plot_points=60, contour=[1,3,5], region=lambda x,y,z: x<=t or y>=t or z<=t)
a = animate([frame(t) for t in srange(.01,1.5,.2)])
a[0]       # long time

a.show()   #

3D rendering not yet implemented
x = var('x')
n = 2
r0 = var('r0')
assume(r0>0)
R = var('R')
assume(R>0)
A = var('A')
assume(A>0)

f(x) = sin(x)^n

v(r0) = integrate(f(x), (x, 0, r0))
V(r0) = R^(n+1) * v(r0) * A
A(r0) = R^n * f(r0) * A

r(x) = solve(x == V(r0), r0)[0].rhs()
I(x) = (A.subs(r0 = r(x))).simplify_full()

show(V)
show(A)
show(r)
show(I)

$\displaystyle r_{0} \ {\mapsto}\ \frac{1}{4} \, A R^{3} {\left(2 \, r_{0} - \sin\left(2 \, r_{0}\right)\right)}$
$\displaystyle r_{0} \ {\mapsto}\ A R^{2} \sin\left(r_{0}\right)^{2}$
$\displaystyle x \ {\mapsto}\ \frac{A R^{3} \sin\left(2 \, r_{0}\right) + 4 \, x}{2 \, A R^{3}}$
$\displaystyle x \ {\mapsto}\ A R^{2} \sin\left(\frac{A R^{3} \cos\left(r_{0}\right) \sin\left(r_{0}\right) + 2 \, x}{A R^{3}}\right)^{2}$

Phi = arctan(x)
Phip = Phi.diff()
Phipp = Phip.diff()

ODE = x * Phipp - 2 * Phip *(Phip - 1)

ODE.full_simplify()

0
f = cos(arccot(x))

f.full_simplify()

x/sqrt(x^2 + 1)
from itertools import imap

(x,y) = var ('x, y')
U = [(1, 2), (3, 4)]
V = [(1, 5), (3, 7)]
f (x,y) = x + y

list(imap(f, [u[1] for u in U], [v[1] for v in V]))

[7, 11]

m = 2
n = 1
G = vector ([cos (m*x), sin (n * x)])

parametric_plot (G, (x, 0,pi/2), aspect_ratio=1)
plot(G [0], (0, 2*pi))
plot(G[1], (0, 2*pi))



f(x) = x^2 + 1
plot(f, (x, 0, 1))

z = var('z')
r = var('r')
l = var('l')
w = var('w')

v = function("u")(z)
vp = v.diff(z)
F = v * vp
G = 2*integrate(F, z)
H = integrate(1/v, z)

Fsin = F.substitute_function(u, sin)
Gsin = G.substitute_function(u, sin)
Hsin = H.substitute_function(u, sin).simplify()
Hsin = -(1/2)*ln((1+cos(z))/(1-cos(z)))
Hsininv = arccos((exp(-2*w) - 1)/(exp(-2*w) + 1))

xi(r) = 2 * arctan(r*l)
xip = xi.diff(r)
xipp = xip.diff(r)

nu(r) = arccos((l^(-2) - r^2)/(l^(-2) + r^2))
nup = nu.diff(r)
nupp = nup.diff(r)

show(Hsin(xi(r)).full_simplify())
show(Hsin(nu(r)).full_simplify())

xiode(r) = (xipp + (1/r) * xip - (1/r^2) * Fsin(xi)).full_simplify()
show(ode)
nuode(r) = (nupp + (1/r) * nup - (1/r^2) * Fsin(nu)).full_simplify()
show(ode)

plot(xi.substitute(l=3), (r,0,1), color="red") + plot(nu.substitute(l=3), (r,0,1), color="blue")

xi(0)
nu(0)

xip(0)
#nup(0)

show(xip.full_simplify())
show(nup.full_simplify())


$\displaystyle -\frac{1}{2} \, \log\left(\frac{1}{l^{2} r^{2}}\right)$
$\displaystyle -\frac{1}{2} \, \log\left(\frac{1}{l^{2} r^{2}}\right)$
$\displaystyle r \ {\mapsto}\ 0$
$\displaystyle r \ {\mapsto}\ 0$
0 0 2*l
$\displaystyle \frac{2 \, l}{l^{2} r^{2} + 1}$
$\displaystyle \frac{2 \, l^{2} r}{{\left(l^{4} r^{4} + 2 \, l^{2} r^{2} + 1\right)} \sqrt{\frac{l^{2} r^{2}}{l^{4} r^{4} + 2 \, l^{2} r^{2} + 1}}}$
n=var('n')
C = n-3
F0 = 2*(n-2)
beta = (1-C)/2 + (1/2)*sqrt((1-C)^2 + 4*F0)
beta


-1/2*n + 1/2*sqrt((n - 4)^2 + 8*n - 16) + 2
function?

File: /ext/sage/sage-8.1/src/sage/calculus/var.pyx
Signature : function(s, *args, **kwds)
Docstring :
Create a formal symbolic function with the name *s*.

INPUT:

* "s" - a string, either a single variable name, or a space or
comma separated list of variable names.

* "**kwds" - keyword arguments. Either one of the following two

keywords can be used to customize latex representation of
symbolic functions:

1. latex_name=LaTeX where "LaTeX" is any valid latex
expression. Ex: f = function('f',
latex_name="mathcal{F}") See EXAMPLES for more.

2. print_latex_func=my_latex_print where
"my_latex_print" is any callable function that returns a
valid latex expression. Ex: f = function('f',
print_latex_func=my_latex_print) See EXAMPLES for an
explicit usage.

Note: The new function is both returned and automatically
injected into the global namespace.  If you use this function in
library code, it is better to use
sage.symbolic.function_factory.function, since it won't touch the
global namespace.

EXAMPLES:

We create a formal function called supersin

sage: function('supersin')
supersin

We can immediately use supersin in symbolic expressions:

sage: y, z, A = var('y z A')
sage: supersin(y+z) + A^3
A^3 + supersin(y + z)

We can define other functions in terms of supersin:

sage: g(x,y) = supersin(x)^2 + sin(y/2)
sage: g
(x, y) |--> supersin(x)^2 + sin(1/2*y)
sage: g.diff(y)
(x, y) |--> 1/2*cos(1/2*y)
sage: k = g.diff(x); k
(x, y) |--> 2*supersin(x)*diff(supersin(x), x)

Custom typesetting of symbolic functions in LaTeX, either using
latex_name keyword:

sage: function('riemann', latex_name="mathcal{R}")
riemann
sage: latex(riemann(x))
mathcal{R}(x)

or passing a custom callable function that returns a latex
expression:

sage: mu,nu = var('mu,nu')
sage: def my_latex_print(self, *args): return "psi_{%s}"%(', '.join(map(latex, args)))
sage: function('psi', print_latex_func=my_latex_print)
psi
sage: latex(psi(mu,nu))
psi_{mu, nu}

In Sage 4.0, you must now use the "substitute_function()" method to
replace functions:

sage: k.substitute_function(supersin, sin)
2*cos(x)*sin(x)

u, v = var('u,v')
alpha = 1/2


u, v = var('u, v', domain='real')

f(v) = exp(v)
F(u, v) = [(1-u) * cos(v), (1-u) * sin(v), f(v=v) * u^2]
M = ParametrizedSurface3D(F,[u,v],'M')

frame = M.natural_frame()
eu = frame[1].substitute({u:0})
ev = frame[2].substitute({u:0})

eu
ev

n = M.normal_vector()
n.substitute({u:0})

h = M.second_fundamental_form_coefficients()
hmatrix = matrix([[h[(1,1)], h[(1,2)]], [h[(2,1)], h[(2,2)]]])
hmatrix.substitute({u:0})


(-cos(v), -sin(v), 0) (-sin(v), cos(v), 0) (0, 0, -1) [-2*e^v 0] [ 0 0]
integrate(cos(x) * cos(x)^2, (x, 0, 2*pi))

0
x, y = var('x, y')
f(x,y)=x^2+y^2
G = plot3d(f, (-1,1), (-2,2))
G.show()
G.save("sphere.png", figsize=1)
G.save("spherezoom.png", zoom=2, figsize=10)
G.save??

3D rendering not yet implemented
   File: /ext/sage/sage-8.1/src/sage/plot/plot3d/base.pyx
Source:
def save(self, filename, **kwds):
"""
Save the graphic in a file.

The file type depends on the file extension you give in the
filename. This can be either:

- an image file (of type: PNG, BMP, GIF, PPM, or TIFF) rendered
using Jmol (default) or Tachyon,

- a Sage object file (of type .sobj) that you can load back later
(a pickle),

- a data file (of type: X3D, STL, AMF, PLY) for export and use in
other software.

For data files, the support is only partial. For instance STL and
AMF only works for triangulated surfaces. The prefered format is X3D.

INPUT:

- filename -- string. Where to save the image or object.

- **kwds -- When specifying an image file to be rendered by Tachyon
or Jmol, any of the viewing options accepted by show() are valid as
keyword arguments to this function and they will behave in the same
way. Accepted keywords include: viewer, verbosity,
figsize, aspect_ratio, frame_aspect_ratio, zoom,
frame, and axes. Default values are provided.

EXAMPLES::

sage: f = tmp_filename(ext='.png')
sage: G = sphere()
sage: G.save(f)

We demonstrate using keyword arguments to control the appearance of the
output image::

sage: G.save(f, zoom=2, figsize=[5, 10])

Using Tachyon instead of the default viewer (Jmol) to create the
image::

sage: G.save(f, viewer='tachyon')

Since Tachyon only outputs PNG images, PIL will be used to convert to
alternate formats::

sage: cube().save(tmp_filename(ext='.gif'), viewer='tachyon')

Here is how to save in one of the data formats::

sage: f = tmp_filename(ext='.x3d')
sage: cube().save(f)

"<Shape><Box size='0.5 0.5 0.5'/><Appearance><Material diffuseColor='0.4 0.4 1.0' shininess='1.0' specularColor='0.0 0.0 0.0'/></Appearance></Shape>"
"""
ext = os.path.splitext(filename)[1].lower()
if ext == '' or ext == '.sobj':
SageObject.save(self, filename)
elif ext in ['.bmp', '.png', '.gif', '.ppm', '.tiff', '.tif',
'.jpg', '.jpeg']:
self.save_image(filename, **kwds)
elif filename.endswith('.spt.zip'):
scene = self._rich_repr_jmol(**kwds)
scene.jmol.save(filename)
elif ext == '.x3d':
with open(filename, 'w') as outfile:
outfile.write(self.x3d())
elif ext == '.stl':
with open(filename, 'wb') as outfile:
outfile.write(self.stl_binary())
elif ext == '.amf':
# todo : zip the output file ?
with open(filename, 'w') as outfile:
outfile.write(self.amf_ascii_string())
elif ext == '.ply':
with open(filename, 'w') as outfile:
outfile.write(self.ply_ascii_string())
else:
raise ValueError('filetype {} not supported by save()'.format(ext))


G.show??

   File: /ext/sage/sage-8.1/src/sage/plot/plot3d/base.pyx
Source:
def show(self, **kwds):
"""
Display graphics immediately

This method attempts to display the graphics immediately,
without waiting for the currently running code (if any) to
return to the command line. Be careful, calling it from within
a loop will potentially launch a large number of external
viewer programs.

INPUT:

-  viewer -- string (default: 'jmol'), how to view
the plot

* 'jmol': Interactive 3D viewer using Java

* 'tachyon': Ray tracer generates a static PNG image

* 'canvas3d': Web-based 3D viewer using JavaScript
and a canvas renderer (Sage notebook only)

* 'threejs': Web-based 3D viewer using JavaScript
and a WebGL renderer

-  verbosity -- display information about rendering
the figure

-  figsize -- (default: 5); x or pair [x,y] for
numbers, e.g., [5,5]; controls the size of the output figure. E.g.,
with Tachyon the number of pixels in each direction is 100 times
figsize[0]. This is ignored for the jmol embedded renderer.

-  aspect_ratio -- (default: "automatic") -- aspect
ratio of the coordinate system itself. Give [1,1,1] to make spheres
look round.

-  frame_aspect_ratio -- (default: "automatic")
aspect ratio of frame that contains the 3d scene.

-  zoom -- (default: 1) how zoomed in

-  frame -- (default: True) if True, draw a
bounding frame with labels

-  axes -- (default: False) if True, draw coordinate
axes

-  **kwds -- other options, which make sense for particular
rendering engines

OUTPUT:

This method does not return anything. Use :meth:save if you
want to save the figure as an image.

CHANGING DEFAULTS: Defaults can be uniformly changed by importing a
dictionary and changing it. For example, here we change the default
so images display without a frame instead of with one::

sage: from sage.plot.plot3d.base import SHOW_DEFAULTS
sage: SHOW_DEFAULTS['frame'] = False

This sphere will not have a frame around it::

sage: sphere((0,0,0))
Graphics3d Object

We change the default back::

sage: SHOW_DEFAULTS['frame'] = True

Now this sphere is enclosed in a frame::

sage: sphere((0,0,0))
Graphics3d Object

EXAMPLES: We illustrate use of the aspect_ratio option::

sage: x, y = var('x,y')
sage: p = plot3d(2*sin(x*y), (x, -pi, pi), (y, -pi, pi))
sage: p.show(aspect_ratio=[1,1,1])

This looks flattened, but filled with the plot::

sage: p.show(frame_aspect_ratio=[1,1,1/16])

This looks flattened, but the plot is square and smaller::

sage: p.show(aspect_ratio=[1,1,1], frame_aspect_ratio=[1,1,1/8])

This example shows indirectly that the defaults
from :func:~sage.plot.plot.plot are dealt with properly::

sage: plot(vector([1,2,3]))
Graphics3d Object

We use the 'canvas3d' backend from inside the notebook to get a view of
the plot rendered inline using HTML canvas::

sage: p.show(viewer='canvas3d')
"""
from sage.repl.rich_output import get_display_manager
dm = get_display_manager()
dm.display_immediately(self, **kwds)


s = sphere()
s.save("sphere.png", xmin=-1, xmax=5)
s.save("sphere.png", zoom=1.2)#, figsize=[2,3])

%latex.add_to_preamble
\renewcommand{\ref}[1]{#1}


%latex

\begin{equation}
\label{eq}
x^2
\end{equation}

See equation $$\eqref{eq}$$

latex.engine()

'pdflatex'
a = var('a')
c = var('c')
x = var('x')
y = function('y')(x)
assume(a*c>0)
desolve(diff(y,x) - a*y^2 + c, y, ivar=x)

1/2*log((a*y(x) - sqrt(a*c))/(a*y(x) + sqrt(a*c)))/sqrt(a*c) == _C + x

t = var('t')
p = var('p')
r = function('r')(t)
r0 = 0
t0 = 0

rp = r.diff(t)

assume(p<0)
de = rp == cot(r)^p

de = rp == r^p

soln = desolve(de, r, ivar=t, ics=[t0, r0])

show(soln)

Error in lines 7-7 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "", line 1, in <module> File "/ext/sage/sage-8.4_1804/local/lib/python2.7/site-packages/sage/symbolic/assumptions.py", line 550, in assume x.assume() File "sage/symbolic/expression.pyx", line 1844, in sage.symbolic.expression.Expression.assume (build/cythonized/sage/symbolic/expression.cpp:14570) raise ValueError("Assumption is %s" % str(s._sage_()[0])) ValueError: Assumption is inconsistent


plot(cot(x), (x, 0.1, pi/2))

t = var('t')
n = var('n')

r(t) = arccos((1/sqrt(2)) * exp(t))
rp(t) = r.diff(t)

rp(t) + cot(r(t))
r(0) == pi/4


0 arccos(1/2*sqrt(2)) == 1/4*pi


R = IntegerRing()
M = MatrixSpace(R, 3, 3)
A = M([538, 282, 273, 549, 270, 273, 546, 273, 273])

show(A)
A.det()
factor(A.det())
show(A.smith_form())

$\displaystyle \left(\begin{array}{rrr} 538 & 282 & 273 \\ 549 & 270 & 273 \\ 546 & 273 & 273 \end{array}\right)$
-819 -1 * 3^2 * 7 * 13
($\displaystyle \left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 273 \end{array}\right)$, $\displaystyle \left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)$, $\displaystyle \left(\begin{array}{rrr} 1 & 9 & -1092 \\ 1 & 8 & -1001 \\ -3 & -26 & 3186 \end{array}\right)$)

D = groups.presentation.Cyclic(8)
F.<u,v> = FreeGroup()
H = F / [u^3, v^3, u*v*u^(-1)*v^(-1)]

a = H.gen(0)
b = H.gen(1)

hom = (D.gens(), [(H.gens(), (b, a*b^2))])
G = D.semidirect_product(H, hom)

delta = G.gen(0)
alpha = G.gen(1)
beta = G.gen(2)

e = delta * alpha

e^2 == delta^2 * beta * alpha
e^3 == delta^3 * alpha^2
e^4 == delta^4 * alpha * beta^2
e^5 == delta^5 * beta^2
e^6 == delta^6 * beta
e^7 == delta^7 * alpha^2 * beta^2
e^8 == G[0]


True True True True True True True
S = SymmetricGroup(8)
subgroups = S.conjugacy_classes_subgroups()
list(map(order, subgroups))

[1, 2, 2, 2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 14, 15, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 18, 18, 18, 18, 18, 18, 18, 20, 20, 20, 21, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 24, 30, 30, 30, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 32, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 36, 40, 42, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 48, 56, 60, 60, 60, 60, 60, 64, 64, 64, 64, 64, 64, 64, 72, 72, 72, 72, 72, 72, 72, 72, 72, 72, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 120, 120, 120, 120, 120, 120, 120, 128, 144, 144, 144, 168, 168, 168, 180, 192, 192, 192, 192, 192, 240, 240, 288, 288, 288, 336, 360, 360, 360, 360, 384, 576, 576, 576, 720, 720, 720, 720, 1152, 1344, 1440, 2520, 5040, 20160, 40320]

p=7
H = S.sylow_subgroup(p)
G = S.quotient(H)
G.order()
G.gens()

p=5
H = S.sylow_subgroup(p)
G = S.quotient(H)
G.order()
G.gens()

p=3
H = S.sylow_subgroup(p)
G = S.quotient(H)
G.order()
G.gens()

p=2
H = S.sylow_subgroup(p)
G = S.quotient(H)
G.order()
G.gens()

Error in lines 3-3 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1188, in execute flags=compile_flags) in namespace, locals File "", line 1, in <module> File "/ext/sage/sage-8.7_1804/local/lib/python2.7/site-packages/sage/groups/perm_gps/permgroup.py", line 2883, in quotient Q = self._gap_() / N._gap_() File "sage/structure/element.pyx", line 1647, in sage.structure.element.Element.__div__ (build/cythonized/sage/structure/element.c:12541) return (<Element>left)._div_(right) File "sage/structure/element.pyx", line 1756, in sage.structure.element.Element._div_ (build/cythonized/sage/structure/element.c:13041) return python_op(other) File "/ext/sage/sage-8.7_1804/local/lib/python2.7/site-packages/sage/interfaces/interface.py", line 1587, in _div_ return self._operation("/", right) File "/ext/sage/sage-8.7_1804/local/lib/python2.7/site-packages/sage/interfaces/interface.py", line 1433, in _operation raise TypeError(msg) TypeError: Gap produced error output Error, <N> must be a normal subgroup of <G> executing \$sage929:=\$sage8 / \\$sage927;;


Subgroup generated by [(1,2,3,4,5,6,7)] of (Symmetric group of order 8! as a permutation group)


S = SymmetricGroup(4)
subgroups = S.conjugacy_classes_subgroups()
subgroups
list(map(order, subgroups))

[Subgroup generated by [()] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(1,3)(2,4)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(3,4)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(2,4,3)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(1,3)(2,4), (1,4)(2,3)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(3,4), (1,2)(3,4)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(1,2)(3,4), (1,3,2,4)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(3,4), (2,4,3)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(3,4), (1,3)(2,4), (1,4)(2,3)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(2,4,3), (1,3)(2,4), (1,4)(2,3)] of (Symmetric group of order 4! as a permutation group), Subgroup generated by [(3,4), (2,4,3), (1,3)(2,4), (1,4)(2,3)] of (Symmetric group of order 4! as a permutation group)] [1, 2, 2, 3, 4, 4, 4, 6, 8, 12, 24]