kill(all)$

maxima : comet approach the sun,problem:quadratic form,part1 version 1.0
peter.vlasschaert@gmail.com,12/05/2022

option :

ratprint:false;

ref:elementary linear algebra,prentice hall international editions ,C.H edwards,1988

a) data : comet approach the sun :
-----------------------------------------
given : 2D data : convert data observation in 3d : 3D → 2D

→ given : 5 points conic .(=Ellipse), pg: 111

cometdata:[[1,9.265],[2,10.892],[3,12.331],[4,13.625],[5,14.799]];

cometdata[1][1];

ordergreat(x,y,z);

principles : Curve Fitting
==================

→ equation : line

s1:A*x+B*y+C=0;

given data:(xi,yi),i=1,2 ,line ,need two points.

s2a:subst([x=x1,y=y1],s1);
s2b:subst([x=x2,y=y2],s1);

→ coef. matrix

s3:coefmatrix([s1,s2a,s2b],[A,B,C]);

→ find :{A,B,C}

method : 1,ratcoef
--------------

s4:determinant(s3)=0;

s4a1:A=ratcoef(lhs(s4),x);
s4a2:B=ratcoef(lhs(s4),y);
s4a3:C=subst([x=0,y=0],lhs(s4));

method :2,differentation
--------------

s4b1:A=diff(lhs(s4),x);
s4b2:B=diff(lhs(s4),y);
s4b3:C=subst([x=0,y=0],lhs(s4));

check : method 1 ? method 2

s5:map(is,[s4a1=s4b1,s4a2=s4b2,s4a3=s4b3]);

b) curve fitting: comet trajectory
--------------------------------------------

→ fitting equation, last position  : 1 , "see above"

t1:A11*x^2+B11*x*y+C11*y^2+D11*x+E11*y+F11=0;

t1a:expand(t1/F11);

t1b:A1*x^2+B1*x*y+C1*y^2+D1*x+E1*y+1=0;

→ coefmatrix

t2b:sort(listofvars(t1b));

t2b1:makelist(ratcoef(lhs(t1b),t2b[i]),i,1,5);

t2b1a:makelist(ratcoef(lhs(t1b)=lhs(t1a),t2b1[i]),i,1,3);

t2b2:makelist(subst([x=cometdata[i][1],y=cometdata[i][2]],t2b1),i,1,5);

t2b3:append([t2b1],t2b2);

t2b4:apply(matrix,t2b3);

→ final :coef matrix "curve fitting"

numer:true$

fpprintprec : 5;

add column to : t2b4

t2b5:addcol(t2b4,[1,1,1,1,1,1]);

see:ref  pg 112

→ find equation: quadratic form

t3c1:determinant(t2b5)=0;

corrected : comet problem , real values *200000

t3c2:expand(t3c1*200000.00);

rem : ≠ ref : pg 112.

c) the implicit plot : quadratic form "t3c2"
------------------------------------------------------

use :plot2d

→ implicit plot

wxplot2d(t3c2, [x, -5, 30], [y, -5, 30],[same_xy, true],[style,[lines]],
 [gnuplot_preamble, "set zeroaxis"]);

see: info , example manual
-----------------------------------

plot2d([[discrete,xy], 2*%pi*sqrt(l/980)], [l,0,50],
        [style, [points,5,2,6], [lines,1,1]],
        [legend,"experiment","theory"],
        [xlabel,"pendulum's length (cm)"], [ylabel,"period (s)"])$

→ implicit plot +comet data

wxplot2d([[discrete,cometdata],t3c2], [x, -5, 30], [y, -5, 30],[same_xy, true],[style,[points,5,2,6],[lines]],[color,red],
 [gnuplot_preamble, "set zeroaxis"],grid2d);

use: wxdraw2d

ref :plot
http://www.scotchildress.com/wxmaxima/Plotting/Plotting.html

→ sun : include point (0,0) = cometdata1

rem : add data to list (a point for the Sun).

cometdata1:append([[0,0]],cometdata);

 wxdraw2d(
   title="comet1-data",
   xaxis=true,
   yaxis=true,
   xaxis_type=solid,
   yaxis_type=solid,
     proportional_axes=xy,
    color =red,
    line_width=1,
   implicit(t3c2,x,-5,25,y,-5,25),
    color = blue,
    point_type=filled_circle,   
   points(cometdata1) 
     
); 

orginal data : cometdata1
---------------------------------

5 points : to find :A1*x^2+B1*y*x+D1*x+C1*y^2+E1*y+1=0,{A1,B1,D1,C1,E1}
1 point  : origin ,sun 

see fig : ref : pg 112, fig 2.9

rem : need to find equations major-axes & minor-axes : call semiaxes "ref: book"

analysis:

rem : vec(v) = vector(v), transpose of vector = transpose(vec(v)) =vT
given :1991.6*x^2-2114.1*y*x-1.8676*10^4*x+1412.0*y^2-4388.1*y-4.4279*10^4=0.0

Q(x,y) = transpose(vec(v))*A*vec(v) =vT*A*v

simp:false;

use :t2b1

z1a:v=[x,y];

z1b:vT= transpose(rhs(z1a));

z1c:A=genmatrix(a, 2, 2);

z1d:Q(x,y)=rhs(z1a).rhs(z1c).rhs(z1b);

simp:true;

z1e:expand(z1d);

A = symmetric matrix (a[i,j])    → A = transpose(A)

z1f:a[2,1]=a[1,2];

z1e1:subst([z1f],z1e);

use :t3c2

z2a:t3c2;

a[1,1]:1991.6;
a[1,2]:-2114.1/2;
a[2,1]:-2114.1/2;
a[2,2]:1412.0;

A: values

z2b:ev(z1c,nouns);

z2b1:determinant(rhs(z2b));

simp:false;

z2c:Q(x,y)=rhs(z1a).rhs(z2b).rhs(z1b);

simp:true;

z2c1:expand(z2c);

rem : xy rotation of the conic (=Ellipse)

Principal Axes : Quadratic Form
---------------------------------------
Q(x',y') =λ[1]*x'^2+λ[2]*y'^2
see : ref ,pg 358

→ find : eigen values ,λ[i],i=1,2

load(eigen)$

z3:eigenvalues(rhs(z2b));

draw axes : minor ,major.

z3a:eigenvalues(rhs(z2b))/25;

→ sqrt

z3b:sqrt(z3a);

z3c:[λ[1]=z3b[1][1],λ[2]=z3b[1][2]];

λ[1] = minoraxes
λ[2] = majoraxes

find : eigenvectors

z4a1:eigenvectors(rhs(z2b));

z4a2:eivects(rhs(z2b));

---------------------------------------------------------------------------------------