#% FILE antonMO.mpl
interface(echo=2);

# based on 
# Example 2 - AR6 Section 8.1 p370; AR7 Section 9.1 p445
 with(linalg,matrix,exponential,vector):
 A:= matrix([ [1, 1],
               [4,-2] ]):
 evalm( exponential(t*A) &* vector([c1-c2/4,c1+c2]) );


# AR6 S8.3; AR7 S9.3 Least squares######################################
# Maple has a function in the linalg package     leastsqrs
# (Get help on it.)
# leastsqrs(M,b)   M matrix   ,  b vector
# finds least squares solution of Mx = b

# Example 1 - AR6 Section 8.3 p391; AR7 Section 9.3 p463
# Find the equation  y = a + b x  for the line of best fit for the points
#   (0,1), (1,3), (2,4), (3,4)
with(linalg,leastsqrs):
M:=array([[1,0],[1,1],[1,2],[1,3]]):
y:=array([1,3,4,4]):
xStar:=leastsqrs(M,y);
# So the line of best fit is  y =  3/2 x + 1
#
# The same problem can be done with functions from the stats package.
# This has some advantages in connection with built-in graphics.
interface(plotdevice=postscript,plotoutput=MOplotfileLsq1);
with(stats,addrecord,statplot,regression):
dat := array([[y,x],[1,0],[3,1],[4,2],[4,3]]):
v:=regression(dat,y=a+b*x);
#                   v := {a = 1.500000000, b = 1.0000000000}
assign(v); # This assigns the above values to a and b
statplot(dat,y=a+b*x);
#
# Other sorts of least-squares fits are possible with the
# regression function.
# (Apologies for the obsolete units used in the text.
# Think of it as a piece of history: "Galileo's data set".)
# Example 3 - AR6 S8.3; AR7 S9.3. With time scaled by factor of 10.
dat := array([[y,x],[-0.18,1],[0.31,2],[1.03,3],[2.48,4],[3.73,5]]):
v:=regression(dat,y=a0+a1*x+a2*x^2);
#  v := {a0 = -.3980000000, a1 = .03471428572, a2 = .1607142857}
#  Run the above for answers, below for a plot.
interface(plotoutput=MOplotfileLsq2);
assign(v); statplot(dat,y=a0+a1*x+a2*x^2); 

# In Maple Vr3
# with(stats): with(fit):
# dat := [[0,1,2,3],[1,3,4,4]]:
# fit[leastsquare[[x,y],y=a+b*x, {a,b}]](dat);
# Giving       {a = 3/2, b = 1}  
#
# Again, other sorts of least-squares fits are possible 
# Example 3 - AR6 S8.3; AR7 S9.3. With time scaled by factor of 10.
# dat := [[1,2,3,4,5],[-0.18,0.31,1.03,2.48,3.73]]:
# fit[leastsquare[[x,y],y=a0+a1*x+a2*x^2 ]](dat);
# Giving as before {a0 = -.3980000000, a1 = .03471428572, a2 =
# .1607142857}


# This approach is to minimize the sum of squares of
# the residuals, using calculus. 
# The quantities for which we solve are {a,b}. 
# Maple Vr1 and Vr2 solve commands do not like being asked
# solve({fa=0,fb=0},{a,b}) 
# but it is relatively straightforward to tell Maple that it is really 
# a simple set of linear equations to solve for a and b, 
# and to use linsolve.
  a:='a': b:='b': f:=sum((y[j]-a-b*x[j])^2,j=1..n):
  fa:=diff(f,a): fb:=diff(f,b):
  with(linalg,hessian,scalarmul,linsolve): 
  A:=hessian(f,[a,b]);
  fa0:=simplify(subs(a=0,b=0,fa));
  fb0:=simplify(subs(a=0,b=0,fb));
  rhvec:=array([-fa0,-fb0]): linsolve(A,rhvec);


# Continuing the least squares example
 with(linalg,transpose,inverse,dotprod,col):
 Mtr:= transpose(M):      evalm( inverse(Mtr&*M) &* Mtr &* y - xStar );
 p:= evalm( M&*xStar):  # This is the projection of  y  into  colSpace(M).
 dotprod( y-p, col(M,1) );                 dotprod( y-p, col(M,2) );

interface(plotoutput=MOplotfile1);
# Here are some graphs of some conic sections
#
# The following is the circle    x^2 + y^2 = 1
 plot( [cos(t), sin(t), t=-Pi..Pi] );

interface(plotoutput=MOplotfile2);
# The following is the ellipse    9*x^2 + 4*y^2 = 36
# (which occurs in Example 4 - AR6 Section 8.6; AR7 Section 9.6)
 plot( [2*cos(t), 3*sin(t), t=-Pi..Pi] );

interface(plotoutput=MOplotfile3);
# The following is the right hand-side of the hyperbola 9*x^2 - 4*y^2 = 36
 plot( [2*sec(t), 3*tan(t), t=-Pi/4..Pi/4] );
# and an alternative parametrisation is as follows
interface(plotoutput=MOplotfile4);
 plot(  [2*cosh(t),3*sinh(t), t=-1 .. 1 ] );

# A more familiar example is the following (rotated) hyperbola
interface(plotoutput=MOplotfile5);
 plot( [x, 1/x, x=-5 .. 5] );


# Moving the origin to the centre of the quadratic locus  Q2  where
 Q2:= a11*X1^2 + a22*X2^2 + 2*a12*X1*X2 + 2*b1*X1 + 2*b2*X2 - c:
 q2:= expand( subs(X1=x1-h,X2=x2-k, Q2) ):
#
# We want the coefficients of x1,x2 to be 0: collect q2 as a polynomial in x1
 qx1:= collect( q2,x1 ):                        # Also as a polynomial in x2:
 qx2:= collect( q2,x2 ):
 coeff1:= coeff(qx1,x1);
 coeff2:= coeff(qx2,x2);
# These include the unwanted term  2 a12 x1 x2:  extract the constant term
 lprint();   solve( {coeff(coeff1,x2,0)=0,coeff(coeff2,x1,0)=0}, {h,k} );
# CHECK that these values give a quadratic in  x1,x2  where the coefficients
#    of x1,x2 are 0.    (This is not as easy as it seems!)


# Example 5 - AR6 Section 8.5 p407; AR7 Section 9.5 p479
 with(linalg):
# a) The symmetric matrix for the quadratic form
#    q(x,y,z)  =  4 x^2 + 4 x y + 4 y^2 + 4 z^2 + 4 x z + 4 y z   is A where
 A:= matrix([ [ 4, 2, 2],
               [ 2, 4, 2],
               [ 2, 2, 4] ]):
# This can be checked as follows
 X:= vector( [x,y,z] ):                    evalm( transpose(X) &* A &* X );
# b) Diagonalize A   
# Example 3 - AR6 Section 8.7 p426; AR7 Section 9.7 p498
 evals:= op(sort([eigenvals(A)]));
  evects:= eigenvects(A);
  if (nops(evects[1][3])=1) then j:=2 else j:=1 fi:
  v1:= evalm(evects[j][3][1]):
  u1:= evalm(v1/norm(v1,2)):
  v2:= evalm(evects[j][3][2]):
  u2:= GramSchmidt([v1,v2])[2]:
  u2:= evalm(u2/norm(u2,2)):
  if (j=1) then
    u3:= evalm(evects[2][3][1]/norm(evects[2][3][1],2)):
  else
    u3:= evalm(evects[1][3][1]/norm(evects[1][3][1],2)):
  fi:
  S:= augment( u1,u2,u3 );
  equal( evalm( transpose(S)&*A&*S ), diag(evals) );
# Thus if the surface q(x,y,z) = C  is rotated so that the x'-, y'-, z'-axes
#  point along  u1, u2, u3  respectively, it becomes the ellipsoid
#
#                2 x'^2 + 2 y'^2 + 8 z'^2 = C .


with(linalg,grad,hessian,eigenvals):
x:='x': y:='y': f:=x^3 + x^4 - x^2*y^2 + (y^4)/2;
g:=grad(f,[x,y]);
cpts:=[solve({g[1]=0,g[2]=0},{x,y})];
h:=hessian(f,[x,y]);
for j from 1 to nops(cpts) do
 print(` ***************************************** `);
 h.j:=subs(cpts[j][1],cpts[j][2],op(h));
 eigenvals(h.j);
 f.j:=subs(cpts[j][1],cpts[j][2],f);
 od;


f:=x^2*y^2  - 2*x -2*y;


f:=2*x*y*z-4*x*z-2*y*z+x^2+y^2+z^2-2*x-4*y+4*z;


#
# Show that matrix Q below is orthog
# and ask Maple to multiply transpose(Q) by Q.
with(linalg,orthog,transpose,multiply): s2:=1/sqrt(2): s3:=1/sqrt(3):
Q:=array([[s3,s2,s2*s3],[s3,0,-2*s2*s3],[s3,-s2,s2*s3]]):
orthog(Q);
#                                   true
multiply(transpose(Q),Q);
# gives the identity matrix verifying that transpose is inverse,
# as characterises real orthog matrices.


