# HW2-Shaowei Sun

## 1916 days ago by bigdata2016

POWER METHOD

global n from numpy import argmax,argmin @interact def f(n=slider(1,20, 1, 7, label="Size of Matrices")): A=random_matrix(ZZ,n) html('<p>This is a %s by %s matrix $A = %s$ created by Sage randomly.</p>'%(n,n,latex(A))) @interact def _(k = input_box(default=QQ(10), label='$k=$'), auto_update=False): x0=vector([1 for i in range(n)]) # Initial guess of eigenvector maxit=20 # Maximum number of iterates dig=8 # number of decimal places to be shown is dig-1 tol=0.0001 err=1 # Initialization of tolerance i=0 while(i<=k): y0=A*x0 ymod=y0.apply_map(abs) imax=argmax(ymod) c1=y0[imax] x1=y0/c1 err=norm(x0-x1) i=i+1 x0=x1 print "Iteration Number:", i-1 print "y"+str(i-1)+"=",y0.n(digits=dig), "c"+str(i-1)+"=", c1.n(digits=dig), "x"+str(i)+"=",x0.n(digits=dig) print "" html('Dominant eigenvalue :$\lambda_{1} \\approx %s$'%(latex(c1.n(digits=dig))))

Size of Matrices

# Newton's Method (From Calculus)

var('n') @interact def graph(g=input_box('sin(x)+1', label="$f(x)=$"), x_range=range_slider(-50, 50, 1, (-50, 50), label='Range of $x$'), size=slider(1, 50, 1, 3, label="Canvas Size"), axis=('Axis', True), initial=input_box('1', label="$x_0 =$"), k=slider(1, 20, 1, 1, label="$k$ times")): f(x)=g # change this to whatever function you like df(x)=diff(f,x); # Sage will compute the derivative of f xn=float(initial) html('          <font color="blue">$f(x)=%s$</font>;    <font color="red">$df(x)=%s$</font><br><br>'%(latex(g),latex(df(x)))) html('$n$    $x_n$') print "====================================================" print 0, " ", xn xm, xma = x_range if axis: axis_TF=1 else: axis_TF=0 p=plot(f, (x, xm, xma) ,color='blue', axes=axis_TF, plot_points=3000, figsize=7.2) NewtonIt(x)=x-(f/df)(x) for i in range(k): p+=line([(xn, 0), (xn, f(xn))], color="red", linestyle='--') p+=plot(df(xn)*(x-xn)+f(xn), (x, xm, xma), color="deeppink", linestyle='-') xn=N(NewtonIt(xn),digits=20) print i+1, " ", xn print "====================================================" show(p)

$f(x)=$
Range of $x$
Canvas Size
Axis
$x_0 =$
$k$ times

# Linear Transformation (From Linear algebra)

var('a, b, c, d, e, f, x, y, x1, y1,k') html('<h3>$f(x,y)=ax^{2}+2bxy+cy^{2}+dx+ey+f=0$<h3>') @interact def qf(a=slider(-50, 50, 1, 34, label="$a$"), b=slider(-50, 50, 1, -12, label="$b$"), c=slider(-50, 50, 1, 41, label="$c$"), d=slider(-50, 50, 1, -40, label="$d$"), e=slider(-50, 50, 1, 20, label="$e$"), f=slider(-50, 50, 1, -25, label="$f$"),size=slider(1, 50, 1, 2, label="Canvas Size"),anima=('Show Process',False)): g(x,y)=a*x^2+2*b*x*y+c*y^2+d*x+e*y+f html('$f(x,y)=%s=0$'%latex(g(x,y))) A=matrix(QQ,2,2,[a,b,b,c]) B=matrix(1,2,[d,e]) xy=matrix(1,2,[x,y]) yx=matrix(2,1,[x,y]) xy1=matrix(1,2,[x1,y1]) yx1=matrix(2,1,[x1,y1]) print ' ' html('<h7><font color="green">Matrix Form: $%s%s%s+%s%s=%s$</font><br><h7>'%(latex(xy),latex(A),latex(yx),latex(B),latex(yx),-f)) print ' ' if det(A)==0: if b==0 and a==0 and c!=0 and d!=0: ychangex=f/d-e^2/(4*c*d) ychangey=e/(2*c) ychange=matrix(2,1,[ychangex,ychangey]) html('$%s=%s+%s$'%(latex(yx1),latex(yx),latex(ychange))) print ' ' html('<h7><font color="red">Standard Form:$y_{1}^2=%s*x_{1}$</font><br><h7> '%latex(-d/c)) p1=implicit_plot(g(x,y),(x,-size, size), (y,-size, size),axes=1,color='red')+implicit_plot(y^2==-d/c*x,(x,-size, size), (y,-size, size),axes=1,color='black',figsize=6.8) if anima==False: show(p1) if anima==True: vv=[] for i in srange(1,11,1): p2=p1+implicit_plot(g(x-ychangex/10*i,y),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(x==-ychangex+ychangex/10*i,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(y==-ychangey,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-ychangex,y),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==-ychangey,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,10,1): p2=p1+implicit_plot(g(x-ychangex,y-ychangey/10*i),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==-ychangey+ychangey/10*i,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) vv.append(p2) for i in srange(1,6,1): p2=p1 vv.append(p2) nh=animate(vv, ymax=size, ymin=-size, xmax=size, xmin=-size,figsize=[5,5]) show(nh) else: if(b==0 and c==0 and a!=0 and e!=0): xchangex=d/(2*a) xchangey=f/e-d^2/(4*a*e) xchange=matrix(2,1,[xchangex,xchangey]) html('$%s=%s+%s$'%(latex(yx1),latex(yx),latex(xchange))) print ' ' html('<h7><font color="red">Standard Form: $x_{1}^2=%s*y_{1}$</font><br><h7> '%latex(-e/a)) p1=implicit_plot(g(x,y),(x,-size, size), (y,-size, size),axes=1,color='red')+implicit_plot(x^2==-e/a*y,(x,-size, size), (y,-size, size),axes=1,color='black') vv=[] if anima==False: show(p1) if anima==True: vv=[] for i in srange(1,11,1): p2=p1+implicit_plot(g(x-xchangex/10*i,y),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(x==-xchangex+xchangex/10*i,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(y==-xchangey,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-xchangex,y),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==-xchangey,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,10,1): p2=p1+implicit_plot(g(x-xchangex,y-xchangey/10*i),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==-xchangey+xchangey/10*i,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,6,1): p2=p1 vv.append(p2) nh=animate(vv, ymax=size, ymin=-size, xmax=size, xmin=-size,figsize=[5,5]) show(nh) else: if (a!=0 and b!=0 and c!=0): va=A.eigenvalues() vavv=A.eigenvectors_right() P1=matrix(2,2,[vavv[0][1][0][0]/vavv[0][1][0].norm(),vavv[1][1][0][0]/vavv[1][1][0].norm(),vavv[0][1][0][1]/vavv[0][1][0].norm(),vavv[1][1][0][1]/vavv[1][1][0].norm()]) P2=matrix(2,2,[vavv[1][1][0][0]/vavv[1][1][0].norm(),vavv[0][1][0][0]/vavv[0][1][0].norm(),vavv[1][1][0][1]/vavv[1][1][0].norm(),vavv[0][1][0][1]/vavv[0][1][0].norm()]) C1=B*P1 C2=B*P2 TF=1 if (va[0]>0): if C1[0][1]!=0: P=P1 vx=va[0] else: if C2[0][0]!=0: P=P2 vy=va[0] TF=2 else: TF=0 html('<h5><font color="red">Error</font><br><h5>') vx2=va[1] vx1=va[0] if (va[0]==0): if C2[0][1]!=0: P=P2 vx=va[1] else: if C1[0][0]!=0: P=P1 vy=va[0] TF=2 else: TF=0 html('<h5><font color="red">Error</font><br><h5>') C=B*P if TF==1: xychange=matrix(2,1,[C[0][0]/(2*vx),f/C[0][1]-(C[0][0])^2/(4*vx*C[0][1])]) html('$%s=%s%s+%s$'%(latex(yx1),latex(P.inverse()),latex(yx),latex(xychange))) print ' ' html('<h7><font color="red">Standard Form:$x_{1}^2=%sy_{1}$</font><br><h7>'%latex(-C[0][1]/vx)) vv=[] dvx=C[0][0]/(2*vx) dvy=f/C[0][1]-(C[0][0])^2/(4*vx*C[0][1]) p1=implicit_plot(g(x,y),(x,-size, size), (y,-size, size),axes=1,color='red')+implicit_plot(vx*(x)^2-C[0][1]*y==0,(x,-size, size), (y,-size, size),axes=1,color='black',figsize=6.8) if P[0][0]==0: k=0 angle=0 else: k=P[1][0]/P[0][0] angle=arctan(k) if anima==False: show(p1) if anima==True: v=[] for i in srange(1,11,1): p2=p1+implicit_plot(g(x-dvx/10*i*cos(angle),y-k*dvx/10*i*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*(y+k*dvx*cot(angle)-k*dvx/10*i*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-dvx*cos(angle),y-k*dvx*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,11,1): p2=p1+implicit_plot(g(x-dvx*cos(angle)-k*dvy/10*cos(angle)*i,y-k*dvx*cos(angle)+dvy/10*cos(angle)*i),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)-k*dvy/10*i*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-dvx*cos(angle)-k*dvy*cos(angle),y-k*dvx*cos(angle)+dvy*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*x,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,10,1): gg(x,y)=g(x-dvx*cos(angle)-k*dvy*cos(angle),y-k*dvx*cos(angle)+dvy*cos(angle)) p2=p1+implicit_plot(gg(x*cos(angle/10*i)-y*sin(angle/10*i),(x)*sin(angle/10*i)+(y)*cos(angle/10*i)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==tan(angle-angle/10*i)*x,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-tan(angle-angle/10*i)*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,6,1): p2=p1 v.append(p2) nh=animate(v, ymax=size, ymin=-size, xmax=size, xmin=-size,figsize=[5,5]) show(nh) if TF==2: yxchange=matrix(2,1,[f/C[0][0]-(C[0][1])^2/(4*vy*C[0][0]),C[0][1]/(2*vy)]) html('$%s=%s%s+%s$'%(latex(yx1),latex(P.inverse()),latex(yx),latex(yxchange))) print ' ' html('<h7><font color="red">Standard Form: $y_{1}^2=%sx_{1}$</font><br><h7>'%latex(-C[0][0]/vy)) vv=[] dvx=C[0][1]/(2*vy) dvy=f/C[0][0]-(C[0][1])^2/(4*vy*C[0][0]) p1=implicit_plot(g(x,y),(x,-size, size), (y,-size, size),axes=1,color='red')+implicit_plot(vy*(y)^2-C[0][0]*x==0,(x,-size, size), (y,-size, size),axes=1,color='black',figsize=6.8) if P[0][0]==0: k=0 angle=0 else: k=P[1][0]/P[0][0] angle=arctan(k) if anima==False: show(p1) if anima==True: for i in srange(1,11,1): p2=p1+implicit_plot(g(x-dvx/10*i*cos(angle),y-k*dvx/10*i*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*(y+k*dvx*cot(angle)-k*dvx/10*i*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-dvx*cos(angle),y-k*dvx*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,11,1): p2=p1+implicit_plot(g(x-dvx*cos(angle)-k*dvy/10*cos(angle)*i,y-k*dvx*cos(angle)+dvy/10*cos(angle)*i),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)-k*dvy/10*i*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-dvx*cos(angle)-k*dvy*cos(angle),y-k*dvx*cos(angle)+dvy*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*x,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,10,1): gg(x,y)=g(x-dvx*cos(angle)-k*dvy*cos(angle),y-k*dvx*cos(angle)+dvy*cos(angle)) p2=p1+implicit_plot(gg(x*cos(angle/10*i)-y*sin(angle/10*i),(x)*sin(angle/10*i)+(y)*cos(angle/10*i)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==tan(angle-angle/10*i)*x,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-tan(angle-angle/10*i)*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') vv.append(p2) for i in srange(1,6,1): p2=p1 vv.append(p2) nh=animate(vv, ymax=size, ymin=-size, xmax=size, xmin=-size,figsize=[5,5]) show(nh) else: html('<h5><font color="red">Error</font><br><h5>') else: va=A.eigenvalues() vavv=A.eigenvectors_right() if va[0]==va[1]: P=matrix(2,2,[1,0,0,1]) else: if a>c: P=matrix(2,2,[vavv[0][1][0][0]/vavv[0][1][0].norm(),vavv[1][1][0][0]/vavv[1][1][0].norm(),vavv[0][1][0][1]/vavv[0][1][0].norm(),vavv[1][1][0][1]/vavv[1][1][0].norm()]) vx=va[0] vy=va[1] else: P=matrix(2,2,[vavv[1][1][0][0]/vavv[1][1][0].norm(),vavv[0][1][0][0]/vavv[0][1][0].norm(),vavv[1][1][0][1]/vavv[1][1][0].norm(),vavv[0][1][0][1]/vavv[0][1][0].norm()]) vx=va[1] vy=va[0] C=B*P D=matrix(2,1,[C[0][0]/(2*vx),C[0][1]/(2*vy)]) f1=f-(C[0][0])^2/(4*vx)-(C[0][1])^2/(4*vy) if f1==0 or (vx>0 and vy>0 and f1>0) or (vx<0 and vy<0 and f1<0): html('<h5><font color="green">Error</font><br><h5>') else: html('$%s=%s%s+%s$'%(latex(yx1),latex(P.inverse()),latex(yx),latex(D))) print ' ' html('<h7><font color="red">Standard Form: $%sx_{1}^2+%sy_{1}^2=%s$</font><br><h7> '%(latex(vx),latex(vy),latex(-f1))) v=[] dvx=C[0][0]/(2*vx) dvy=C[0][1]/(2*vy) p1=implicit_plot(g(x,y),(x,-size, size), (y,-size, size),axes=1,color='red')+implicit_plot(vx*(x)^2+vy*(y)^2+f1==0,(x,-size, size), (y,-size, size),axes=1,color='black',figsize=6.8) if P[0][0]==0: k=0 angle=0 else: k=P[1][0]/P[0][0] angle=arctan(k) if anima==False: show(p1) if anima==True: for i in srange(1,11,1): p2=p1+implicit_plot(g(x-dvx/10*i*cos(angle),y-k*dvx/10*i*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*(y+k*dvx*cot(angle)-k*dvx/10*i*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-dvx*cos(angle),y-k*dvx*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,11,1): p2=p1+implicit_plot(g(x-dvx*cos(angle)-k*dvy/10*cos(angle)*i,y-k*dvx*cos(angle)+dvy/10*cos(angle)*i),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*(x+k*dvy*cot(angle)-k*dvy/10*i*cot(angle)),(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,6,1): p2=p1+implicit_plot(g(x-dvx*cos(angle)-k*dvy*cos(angle),y-k*dvx*cos(angle)+dvy*cos(angle)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==k*x,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-k*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) v.append(p2) for i in srange(1,10,1): gg(x,y)=g(x-dvx*cos(angle)-k*dvy*cos(angle),y-k*dvx*cos(angle)+dvy*cos(angle)) p2=p1+implicit_plot(gg(x*cos(angle/10*i)-y*sin(angle/10*i),(x)*sin(angle/10*i)+(y)*cos(angle/10*i)),(x,-size, size), (y,-size, size),axes=1) p2=p2+implicit_plot(y==tan(angle-angle/10*i)*x,(x,-size, size), (y,-size, size),axes=1,linestyle='--')+implicit_plot(x==-tan(angle-angle/10*i)*y,(x,-size, size), (y,-size, size),axes=1,linestyle='--') v.append(p2) for i in srange(1,6,1): p2=p1 v.append(p2) nh=animate(v, ymax=size, ymin=-size, xmax=size, xmin=-size,figsize=[5,5]) show(nh)

## Click to the left again to hide and once more to show the dynamic interactive window

Simple Regression Analysis (From Statistics)

%r IQ=c(110, 130, 125, 120, 115, 120, 125, 130, 150, 140, 100, 110, 115, 120, 135) Math=c(75, 90, 80, 80, 70, 75, 90, 95, 90, 85, 60, 65, 75, 75, 90) x=IQ y=Math lsm=lm(y~x) summary(lsm) plot(IQ, Math, col="blue", pch=16, main="Scatter Diagram and Fitted Regression Line ") abline(coef(lsm), col="red") dev.off()
 Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -7.824 -4.243 -1.012 3.213 10.626 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.0495 14.9894 -0.203 0.842 x 0.6725 0.1213 5.546 9.45e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 5.803 on 13 degrees of freedom Multiple R-squared: 0.7029, Adjusted R-squared: 0.6801 F-statistic: 30.76 on 1 and 13 DF, p-value: 9.447e-05 null device 1  Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -7.824 -4.243 -1.012 3.213 10.626 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -3.0495 14.9894 -0.203 0.842 x 0.6725 0.1213 5.546 9.45e-05 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 5.803 on 13 degrees of freedom Multiple R-squared: 0.7029, Adjusted R-squared: 0.6801 F-statistic: 30.76 on 1 and 13 DF, p-value: 9.447e-05 null device 1