Computing with SageMath Directly from the Forum

[Eli]

Type some Sage code below and press Compute (Test examples). Choose a different language from the dropdown menu (say Macaulay2), write some code in the editor, and click the "Compute" button to run it (reload browser if you do not see the editor).

[Eli]

Try to integrate 1/(1+x^4) with Sage or Maxima, i.e., integrate(1/(1 + x^4), x);

The answer should be

\begin{equation}
\begin{split}
& \int\frac{1}{1 + x^4}{\rm d}x \\
& = 1/4\sqrt{2}\times{\rm tan}^{-1}\left(1/2\sqrt{2}\times(2x + \sqrt{2})\right) + 1/4\sqrt{2}\times {\rm tan}^{-1}\left(1/2\times \sqrt{2}\times(2x - \sqrt{2})\right)\\
& + 1/8\sqrt{2}\times{\rm log}(x^2 + \sqrt{2}\times x + 1) - 1/8\sqrt{2}\times{\rm log}(x^2 - \sqrt{2}\times x + 1)
\end{split}
\end{equation}


Try this in Macaulay2:

Code: [Select all] [Expand/Collapse]

1. R = QQ[x,y,z]
2. curve = ideal( x^4-y^5, x^3-y^7 )
3. gb curve
4. dim curve
5. degree curve
6. curve1 = saturate(curve,ideal(x))

[Eli]

Test these solutions in Sage (see reference manual here):

Code: [Select all] [Expand/Collapse]

1. x,y,z=var('x,y,z')
2. # Next we define the parameters
3. sigma=10
4. rho=28
5. beta=8/3
6. # The Lorenz equations
7. lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
8. # Time and initial conditions
9. times=srange(0,50.05,0.005)
10. ics=[0,1,1]
11. sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14)
12. vis=line(zip(sol[:,0],sol[:,1]))
13. vis.show()

The output is:

Code: [Select all] [Expand/Collapse]

1. #from sage.calculus.desolvers import desolve_odeint -- optional
2. x,y=var('x,y')
3. f=[x*(1-y),-y*(1-x)]
4. sol=desolve_odeint(f,[0.5,2],srange(0,10,0.1),[x,y])
5. vis=line(zip(sol[:,0],sol[:,1]))
6. print(sol)
7. vis.show()

Code: [Select all] [Expand/Collapse]

1. #y = function('y')(x)
2. #eq = diff(y,x) - (x - 2)/(y^2 + 1)*sin(x - 2) - sin(y) == 0
3. #desolve_laplace(eq,y,ics=[0,3])
4. y1,y2,y3=var('y1,y2,y3')
5. f1=77.27*(y2+y1*(1-8.375*1e-6*y1-y2))
6. f2=1/77.27*(y3-(1+y1)*y2)
7. f3=0.16*(y1-y3)
8. f=[f1,f2,f3]
9. ci=[0.2,0.4,0.7]
10. t=srange(0,10,0.01)
11. v=[y1,y2,y3]
12. sol=desolve_odeint(f,ci,t,v,rtol=1e-3,atol=1e-4,h0=0.1,hmax=1,hmin=1e-4,mxstep=1000,mxords=17)
13. print(sol)

Solve system of ODEs:

Code: [Select all] [Expand/Collapse]

1. t = var('t')
2. x = function('x')(t)
3. y = function('y')(t)
4. de1 = diff(x,t) + y - 1 == 0
5. de2 = diff(y,t) - x + 1 == 0
6. sol = desolve_system([de1, de2], [x,y])
7. #Add ics if you want
8. #sol = desolve_system([de1, de2], [x,y],ics=[0,0.5,1],ivar=t)
9. print(sol)

Code: [Select all] [Expand/Collapse]

1. from sage.calculus.desolvers import desolve_system_rk4
2. x,y,t=var('x y t')
3. P=desolve_system_rk4([x*(1-y),-y*(1-x)],[x,y],ics=[0,0.5,2],ivar=t,end_points=20)
4. Q=[ [i,j] for i,j,k in P]
5. LP=list_plot(Q)
6.  
7. Q=[ [j,k] for i,j,k in P]
8. LP=list_plot(Q)
9. plot(LP)

Code: [Select all] [Expand/Collapse]

1. var('t,x,y,X,Y')
2. f(t,x,y,X,Y)=[X, Y, -x/(x^2+y^2)^(3/2), -y/(x^2+y^2)^(3/2)]
3. ics = [0.8, 0, 0, 1.22474487139159]
4. t = 100*pi
5. sol = desolve_mintides(f, ics, 0, t, t, 1e-12, 1e-12) # optional -tides
6. sol # optional -tides # abs tol 1e-5

[Eli]

Plotting in Sage (see more from here):

Code: Select all

f(x) = (x-3)*(x-5)*(x-7)+40
P = line([(2,0),(2,f(2))], color='black')
P += line([(8,0),(8,f(8))], color='black')
P += polygon([(2,0),(2,f(2))] + [(x, f(x)) for x in [2,2.1,..,8]] + [(8,0),(2,0)],  rgbcolor=(0.8,0.8,0.8),aspect_ratio='automatic')
P += text("$\\int_{a}^b f(x) dx$", (5, 20), fontsize=16, color='black')
P += plot(f, (1, 8.5), thickness=3)
P    # show the result

Code: [Select all] [Expand/Collapse]

1. var("x y")
2. G = Graphics()
3. counter = 0
4. for col in colors.keys():  # Long time
5.     G += implicit_plot(x^2 + y^2 == 1 + counter*.1, (x,-4,4),(y,-4,4), color=col)
6.     counter += 1
7. G

Code: [Select all] [Expand/Collapse]

1. var("x y")
2. contour_plot(y^2 + 1 - x^3 - x, (x,-pi,pi), (y,-pi,pi), fill=False, cmap='hsv', labels=True)

Code: [Select all] [Expand/Collapse]

1. var("x y")
2. parametric_plot([cos(x) + 2 * cos(x/4), sin(x) - 2 * sin(x/4)], (x,0, 8*pi), fill=True)

Code: [Select all] [Expand/Collapse]

1. def b(n): return lambda x: bessel_J(n, x)
2. plot([b(n) for n in [1..5]], 0, 20, fill='axis')

See the attached outputs.

[Eli]

Try to solve

\begin{equation}
\frac{{\rm d}y}{{\rm d}x} = \frac{x - 2}{y^2 + 1}\times {\rm sin}(x - 2) - {\rm sin}(y).
\end{equation}

[Eli]

We can label and decorate our plots:

Code: [Select all] [Expand/Collapse]

1. x,y,z=var('x,y,z')
2. # Next we define the parameters
3. sigma=10
4. rho=28
5. beta=8/3
6. # The Lorenz equations
7. lorenz=[sigma*(y-x),x*(rho-z)-y,x*y-beta*z]
8. # Time and initial conditions
9. times=srange(0,50.05,0.005)
10. ics=[0,1,1]
11. sol=desolve_odeint(lorenz,ics,times,[x,y,z],rtol=1e-13,atol=1e-14)
12. vis=line(zip(sol[:,0],sol[:,1]), color='purple')
13. vis.axes_labels(['X axis', 'Y axis'])
14. vis.legend(False)
15. vis.show(title='Lorenz Equations', frame=True, legend_loc="lower right")
16. vis.show()

[Eli]

Try the following Algebra problems in Sage:

Example 1: Transpose of a Matrix

Code: [Select all] [Expand/Collapse]

1. A = matrix(QQ, 5, range(25))
2. A.T

Example 2: LU Factorization of a Matrix

Code: [Select all] [Expand/Collapse]

1. A = matrix(QQ, [[ 2, 1, 1],
2.                 [ 4,-6, 0],
3.                 [-2, 7, 2]])
4.  
5. (P,L,U) = A.LU()
6. #table([[A, '=', P, '$\\cdot$', L, '$\\cdot$', U]])
7. print(P, L, U, sep='\n')

Code: [Select all] [Expand/Collapse]

1. A = matrix(QQ, [[ 2, 1, 1],
2.                 [ 4,-6, 0],
3.                 [-2, 7, 2]])
4.  
5. (P,L,U) = A.LU(pivot='nonzero')
6. #table([[A, '=', P, '$\\cdot$', L, '$\\cdot$', U]])
7. print(P, L, U, sep='\n')

Code: [Select all] [Expand/Collapse]

1. A = matrix(QQ, [[ 1, -4,  1,  0, -2,  1, 3,  3,  2],
2. [-1,  4,  0, -4,  0, -4, 5, -7, -7],
3. [ 0,  0,  1, -4, -1, -3, 6, -5, -6],
4. [-2,  8, -1, -4,  2, -4, 1, -8, -7],
5. [ 1, -4,  2, -4, -3,  2, 5,  6,  4]])
6. P, L, U = A.LU()
7. U
8. A.rref()
9. A.pivots()

Code: [Select all] [Expand/Collapse]

1. D = matrix(QQ, [[ 1,  0,  2,  0, -2, -1],
2.  [ 3, -2,  3, -1,  0,  6],
3.  [-4,  2, -3,  1, -1, -8],
4.  [-2,  2, -3,  2,  1,  0],
5.  [ 0, -1, -1,  0,  2,  5],
6. [-1,  2, -4, -1,  5, -3]])
7. P, L, U = D.LU(pivot='nonzero')
8. P
9. L
10. U
11. D == L*U

Example 3: Modified Gram-Schmidt Algorithm -- QR decomposition

Code: [Select all] [Expand/Collapse]

1. A = matrix(QQbar, [[-2, 0, -4, -1, -1],
2.  [-2, 1, -6, -3, -1],
3. [1, 1, 7, 4, 5],
4.  [3, 0, 8, 3, 3],
5.  [-1, 1, -6, -6, 5]])
6. Q, R = A.QR()
7. Q
8. R
9. Q.conjugate_transpose()*Q
10. Q*R == A

See more examples here.

[Eli]

Knots Using Sage:

Code: [Select all] [Expand/Collapse]

1. B = BraidGroup(4)
2. K = Knot(B([1,1,1,2,-1,2,-3,2,-3]))
3. plot(K)
4. K.alexander_polynomial()
5. K.jones_polynomial()
6. K.determinant()
7. K.signature()

[Eli]

Interactive Fourier Series with Sage (see source):

Code: [Select all] [Expand/Collapse]

1. ## Interactive Fourier Series, via: http://www.walkingrandomly.com/?p=1879
2. def ftermSquare(n):
3.  return(1/n*sin(n*x*pi/3))
4.  
5. def ftermSawtooth(n):
6.  return(1/n*sin(n*x*pi/3))
7.  
8. def ftermParabola(n):
9.  return((-1)^n/n^2 * cos(n*x))
10.  
11. def fseriesSquare(n):
12.  return(4/pi*sum(ftermSquare(i) for i in range (1,2*n,2)))
13.  
14. def fseriesSawtooth(n):
15.  return(1/2-1/pi*sum(ftermSawtooth(i) for i in range (1,n)))
16.  
17. def fseriesParabola(n):
18.  return(pi^2/3 + 4*sum(ftermParabola(i) for i in range(1,n)))
19.  
20. @interact
21. def plotFourier(n=slider(1, 30,1,3,'Number of terms')
22. ,Function=['Square Wave','Saw Tooth','Periodic Parabola']):
23.     if Function=='Saw Tooth':
24.      show(plot(fseriesSawtooth(n),x,-6,6,figsize=(7,3)))
25.     if Function=='Square Wave':
26.      show(plot(fseriesSquare(n),x,-6,6,figsize=(7,3)))
27.     if Function=='Periodic Parabola':
28.      show(plot(fseriesParabola(n),x,-6,6,figsize=(7,3)))