CoCalc -- 2603.sagews

All published worksheets from http://sagenb.org

Project: sagenb.org published worksheets

Views: ¹⁶⁸⁷³¹
Image: ubuntu2004

x=0;x+=1;x

1

x=2;x**=2;x

4

x=[123,'xyz'];x+=[011];x

[123, 'xyz', 9]

x=y=z=123

123

123

123

x,y,z=1,2,3

1

2

3

x=666;y=888

x,y=y,x

888

666

n=var('n')

a_n=(1+1/n)**n

a_n.limit(n=oo)

e

limit(sin(n*pi/3),n=oo)

ind

n=var('n')

limit((n**2-1)/(n+5),n=oo)

+Infinity

x=var('x')

diff(sin(x),x)

cos(x)

diff(sin(x),x,2)

-sin(x)

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

f=x^2*y+x*exp(y)

f.diff(x,y)

2*x + e^y

var('x,a')

(x, a)

f=exp(sin(a-x^2))/x

f.derivative(x)

-2*e^(sin(-x^2 + a))*cos(-x^2 + a) - e^(sin(-x^2 + a))/x^2

x=var('x')

f=sin(x)

taylor(f,x,0,6)

1/120*x^5 - 1/6*x^3 + x

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

f=sin(x*y)

taylor(f,(x,0),(y,-1),4)

x=var('x')

-1/2*(y + 1)*x^3 + 1/6*x^3 + (y + 1)*x - x

x=var('x')

integral(x*sin(x**2),x)

-1/2*cos(x^2)

integral(x^2,x,0,1)

1/3

x=var('x')

integral(x**(-2),1,infinity)

1

var('x,k,w')

(x, k, w)

f=x^3*exp(k*x)*sin(w*x)

f.integrate(x)

-(((k^6*w + 3*k^4*w^3 + 3*k^2*w^5 + w^7)*x^3 - 24*k^3*w + 24*k*w^3 - 6*(k^5*w + 2*k^3*w^3 + k*w^5)*x^2 + 6*(3*k^4*w + 2*k^2*w^3 - w^5)*x)*e^(k*x)*cos(w*x) - ((k^7 + 3*k^5*w^2 + 3*k^3*w^4 + k*w^6)*x^3 - 6*k^4 + 36*k^2*w^2 - 6*w^4 - 3*(k^6 + k^4*w^2 - k^2*w^4 - w^6)*x^2 + 6*(k^5 - 2*k^3*w^2 - 3*k*w^4)*x)*e^(k*x)*sin(w*x))/(k^8 + 4*k^6*w^2 + 6*k^4*w^4 + 4*k^2*w^6 + w^8)

m=var('m')

sum(1/(m*(m+1)),m,1,oo)

1

sum(1/(m*(m+1)),m,1,10)

10/11

sum(1/(m*(m+1)),m,1,100)

100/101

t=var('t')

x=function('x',t)

a,b=var('a,b')

DE=a*diff(x,t)+b*x-1

desolve(DE,[x,t])

(c + e^(b*t/a)/b)*e^(-b*t/a)

plot(sin,(0,2*pi),color='red')

P1=PolynomialRing(ZZ,'x')

x=P1.gen()

P=-1*x^5+3*x^8-1+6*x^2;P

3*x^8 - x^5 + 6*x^2 - 1

P2=PolynomialRing(RR,'t')

t=P2.gen()

q=2*t^2+5*t^4-4;q

5.00000000000000*t^4 + 2.00000000000000*t^2 - 4.00000000000000

P(2)

759

y=matrix(RR,[[1,2],[3,4]]);y

[1.00000000000000 2.00000000000000]
[3.00000000000000 4.00000000000000]

P(y)

[  496225.000000000   723212.000000000]
[1.08481800000000e6 1.58104300000000e6]

P.change_ring(QQ)

3*x^8 - x^5 + 6*x^2 - 1

q.change_ring(QQ)

5*t^4 + 2*t^2 - 4

P.roots()

[]

P.complex_roots()

[-0.405551899600305, 0.410122280158637, -0.931044157794464 - 0.617070397270502*I, -0.931044157794464 + 0.617070397270502*I, -0.0409077197884951 - 1.14412829588053*I, -0.0409077197884951 + 1.14412829588053*I, 0.969666687303793 - 0.534135779668537*I, 0.969666687303793 + 0.534135779668537*I]

q1=(x-1)*(x-2)*(x-3)^2*(x-4);q

5.00000000000000*t^4 + 2.00000000000000*t^2 - 4.00000000000000

q1.roots()

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

P.is_irreducible()

True

q2=x^3-9*x^2+20*x-12;q2

x^3 - 9*x^2 + 20*x - 12

gcd(q1,q2)

x^2 - 3*x + 2

lcm(q1,q2)

x^6 - 19*x^5 + 143*x^4 - 545*x^3 + 1104*x^2 - 1116*x + 432

q1.factor()

(x - 4) * (x - 2) * (x - 1) * (x - 3)^2

q2.factor()

(x - 6) * (x - 2) * (x - 1)

r.lm#

lm

x1=[1,1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5,6,6,6,6,7,8,8,8,8,10,10,10,10,11,11,12,12,13,13,14,15,16,16,17,20]

x2=[1,0,1,0,0,1,0,0,0,0,1,1,1,0,1,0,0,0,0,1,0,1,0,1,1,0,1,1,0,0,0,1,1,1,0,0,1,0,1,0,1,1,0,0,0,0]

x3=[1,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,1,1,0,0,0,1,0,1,0,1,1,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,1]

x4=[0,0,0,1,0,1,1,0,0,1,0,1,0,0,0,0,1,1,0,0,0,0,1,1,0,1,0,0,0,0,1,0,1,1,0,1,0,0,1,1,0,1,1,0,1,0]

yy=[13876,11608,18701,11283,11767,20872,11772,10535,12195,12313,14975,21371,19800,11417,20263,13231,12884,13245,13677,15965,12366,21352,13839,22884,16978,14803,17404,22184,13548,14467,15942,23174,23780,25410,14861,16882,24170,25990,26330,17949,25685,27837,18838,17483,19207,19346]

x=r.matrix([x1,x2,x3,x4],46,4)

y=r.matrix(yy)

lm=r.lm(y.tilde(x));lm

Call:
lm(formula = sage142)

Coefficients:
(Intercept)     sage1911     sage1912     sage1913     sage1914  
    12024.5        485.5        418.6       1894.2       2905.3

r.confint(lm)

2.5 %     97.5 %
(Intercept)  8951.7641 15097.2787
sage1911      237.8817   733.1953
sage1912    -2112.8437  2950.0909
sage1913    -1509.6394  5298.1177
sage1914     -135.0331  5945.6063

residuals=r.residuals(lm)

resid_data=residuals.sage()['DATA']

resid_x1_list=[(x1[k],resid_data[k]) for k in range(len(x1))]

list_plot(resid_x1_list)

resid_x2_list=[(x2[k],resid_data[k]) for k in range(len(x1))]

list_plot(resid_x2_list)

x24=[0 for k in range(len(x2))]

for k in range(len(x24)):
   if x2[k]==0:
      if x3[k]==1:
          x24[k]=1
      elif x4[k]==1:
          x24[k]=2
      else:
         x24[k]=3
   else:
       if x3[k]==1:
           x24[k]=4
       elif x4[k]==1:
           x24[k]=5
       else:
          x24[k]=6

resid_x24_list=[(x24[k],resid_data[k]) for k in range(len(x24))]

list_plot(resid_x24_list)

x5=[x2[k]*x3[k] for k in range(len(x2))]

x6=[x2[k]*x4[k] for k in range(len(x2))]

x=r.matrix([x1,x2,x3,x4,x5,x6],46,6)

y=r.matrix(yy)

lm=r.lm(y.tilde(x));lm

Call:
lm(formula = sage235)

Coefficients:
(Intercept)     sage2881     sage2882     sage2883     sage2884     sage2885     sage2886  
    12582.2        483.9       -473.2       1835.8       1766.7       -687.1       2374.9

for k in range(len(resid_data)):
   if resid_data[k]==min(resid_data):
      error_station=k

error_station

43

x1.pop(error_station)

17

x2.pop(error_station)

0

x3.pop(error_station)

1

x4.pop(error_station)

0

x5.pop(error_station)

0

x6.pop(error_station)

0

yy.pop(error_station)

17483

x=r.matrix([x1,x2,x3,x4,x5,x6],45,6)

y=r.matrix(yy)

lm=r.lm(y.tilde(x));lm

Call:
lm(formula = sage245)

Coefficients:
(Intercept)     sage2871     sage2872     sage2873     sage2874     sage2875     sage2876  
    12457.5        531.9       -696.3       1962.0       1676.9       -729.3       2458.8

r.confint(lm)

2.5 %     97.5 %
(Intercept)  8688.6768 16226.2917
sage2871      261.2944   802.4869
sage2872    -5581.0127  4188.3795
sage2873    -3031.5392  6955.4567
sage2874    -2938.4118  6292.1197
sage2875    -7968.5830  6510.0328
sage2876    -3733.5680  8651.1374

residuals=r.residuals(lm)

resid_data=residuals.sage()['DATA']

resid_x24_list=[(x24[k],resid_data[k]) for k in range(len(x24))]

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_99.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("cmVzaWRfeDI0X2xpc3Q9Wyh4MjRba10scmVzaWRfZGF0YVtrXSkgZm9yIGsgaW4gcmFuZ2UobGVuKHgyNCkpXQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
  File "", line 1, in <module>
    
  File "/tmp/tmpzaseTZ/___code___.py", line 2, in <module>
    exec compile(u'resid_x24_list=[(x24[k],resid_data[k]) for k in range(len(x24))]' + '\n', '', 'single')
  File "", line 1, in <module>
    
IndexError: list index out of range

list_plot(resid_x24_list)

for k in range(len(resid_data)):
   if resid_data[k]==min(resid_data):
      error_station=k

error_station

43

x1.pop(error_station)

20

x2.pop(error_station)

0

x3.pop(error_station)

0

x4.pop(error_station)

1

x5.pop(error_station)

0

x6.pop(error_station)

0

yy.pop(error_station)

19207

x=r.matrix([x1,x2,x3,x4,x5,x6],45,6)

y=r.matrix(yy)

lm=r.lm(y.tilde(x));lm

[1] 0

r.confint(lm)

[1] 1

def buffontest(1,t,n):
    var('x')
    p=parametric_plot((0,x),(x,-2,2),color="red",frame=True)
    p+=parametric_plot((0,x),(x,-2,2),color="red",frame=True)
    crossnumber=0
    for k in range(n):
        degree=random()*pi
        width=random()*t
        temp=random()
        x1=width-1/2*sin(degree)
        x2=width+1/2*sin(degree)
        y1=temp-1/2*cos(degree)
        y2=temp+1/2*cos(degree)
        if width>t/2:
           width=t-width
        if width-1/2*sin(degree)<0:
           crossnumber+=1
           p+=line([(x1,y1),(x2,y2)],color="#00FF00",frame=True)
        else:
           p+=line([(x1,y1),(x2,y2)],frame=True)
    p.show(aspect_ratio=1,xmin=-0.5,xmax=2.5,ymin=-0.5,ymax=1.5)
    return 2*1*n/t/crossnumber

def motecarlo_integral(s,x,a,b,m):
    myintegral=0.0
    miny_x=minimize(s(x),[(a+b)/2.0],disp=0)[0]
    miny=s(miny_x)
    maxy_x=minimize(-s(x),[(a+b)/2.0],disp=0)[0]
    miny=s(miny_x)
    p=plot(a(x),(x,a,b))
    for k in range(m):
       rnd1=random()*abs(a-b)+min(a,b)
       rnd2=random()*(maxy-miny)+miny
       y=s(x=rnd1)
       if 0<rnd2<n(y):
          myintegral+=1.0
          p+=scatter_plot([(rnd1,rnd2)],marker='o',facecolor='red')
       elif y.n()<rnd2<0:
          myintegral-=1.0
          p+=scatter_plot([(rnd1,rnd2)],marker='s',facecolor='green')
       else:
          p+=list_plot([(rnd1,rnd2)])
    p.show(aspect_ratio=1)
    return n(myintegral/m*(max(a,b)-min(a,b))*(maxy-miny))