gad.m

这是几何代数的matlab工具包

function GAD(GAn)
%GAD: run sample code.

try
if ( GAn == 1 ) 
     GAps = '>>>> ';
     disp('>> % 	VECTORS');
     % 	VECTORS
     GAfigure; clc; %/
     disp('>> % 	VECTORS');
     % 	VECTORS
     global v w; %/
     clf; %/
     disp('>> %	Vectors can be defined using coordinates:');
     %	Vectors can be defined using coordinates:
     disp('>> v = e1 + e2;');
     v = e1 + e2;
     disp('>> w = e2 + 2*e3;');
     w = e2 + 2*e3;
     disp('>> draw(v,''b'');');
     draw(v,'b');
     GAtext(0.5*v -0.1*unit(w),'v'); %/
     disp('>> draw(w,''g'');');
     draw(w,'g');
     GAtext(0.5*w-0.1*unit(v),'w'); %/ 
     axis('vis3d'); %/
     va = [-0.1 1 0 2 -0.1 2]; %/
     axis(va); %/
     GAprompt; %/
     GAorbiter(60,6); %/
     GAprompt; %/
     disp('>> % 	Adding vectors, as usual:');
     % 	Adding vectors, as usual:
     disp('>> v+w');
     v+w
     draw(v+w,'r'); %/
     GAtext(0.7*(v+w) -0.1*unit(w),'v+w'); %/
     DrawPolyline({v,v+w},'k'); %/
     DrawPolyline({w,v+w},'k'); %/
     disp('>> % 	But vector addition is coordinate-free,');
     % 	But vector addition is coordinate-free,
     disp('>> %	so we will drop the axes.');
     %	so we will drop the axes.
     axis(va); %/
     GAprompt; %/
     axis off; %/
     GAprompt; %/
     disp('>> %	Vectors `are'' geometrical objects in space.');
     %	Vectors `are' geometrical objects in space.
     GAorbiter(360,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 2 ) 
     GAps = '>>>> ';
     disp('>> % 	SPANNING: THE OUTER PRODUCT');
     % 	SPANNING: THE OUTER PRODUCT
     GAfigure; clc; %/
     disp('>> % 	SPANNING: THE OUTER PRODUCT');
     % 	SPANNING: THE OUTER PRODUCT
     disp('>> %');
     %
     global v w; %/
     clf; %/
     disp('>> % 	The outer product is anti-symmetric:');
     % 	The outer product is anti-symmetric:
     disp('>> %');
     %
     fprintf(1,'>> e1^e2  ');
     input('');
     e1^e2 %w
     fprintf(1,'>> e2^e1  ');
     input('');
     e2^e1 %w
     fprintf(1,'>> e1^e1  ');
     input('');
     e1^e1 %w
     GAprompt; %/
     disp('>> % 	and the outer product is linear:');
     % 	and the outer product is linear:
     fprintf(1,'\n');     fprintf(1,'>> e1^(e2+e3)  ');
     input('');
     e1^(e2+e3) %w
     disp('>> %	All bivectors expressible on BIVECTOR BASIS.');
     %	All bivectors expressible on BIVECTOR BASIS.
     fprintf(1,'\n');     disp('>> % 	Those properties yield certain identities:');
     % 	Those properties yield certain identities:
     fprintf(1,'\n');     fprintf(1,'>> e1^(e1+e2)  ');
     input('');
     e1^(e1+e2) %w
     GAprompt; %/
     disp('>> %');
     %
     disp('>> % Now let us draw these bivectors');
     % Now let us draw these bivectors
     fprintf(1,'\n');     disp('>> v = e1 + e2;');
     v = e1 + e2;
     disp('>> w = e2 + 2*e3;');
     w = e2 + 2*e3;
     disp('>> B = v^w;');
     B = v^w;
     disp('>> DrawBivector(v,w);');
     DrawBivector(v,w);
     GAview([30 30]); %/
     GAtext(0.5*v-0.1*unit(w)+0.1*unit(grade(((v^w)/I3),1)),'v'); %/
     GAtext(0.5*w+0.1*unit(v)+0.1*unit(grade(((v^w)/I3),1))+v,'w'); %/
     GAtext(0.5*v+0.7*w-0.2*unit((v^w)*I3),'v \wedge w'); %/
     axis('vis3d'); %/
     GAprompt; %/ 
     disp('>> %	Again, the actual object is coordinate-free.');
     %	Again, the actual object is coordinate-free.
     axis off; %/
     GAprompt; %/
     GAorbiter(380,10); %/
     disp('>> % 	The bivector v^w has dimension, direction, ');
     % 	The bivector v^w has dimension, direction, 
     disp('>> %	sense and area, but NO SHAPE');
     %	sense and area, but NO SHAPE
     GAprompt; %/
     va = [-1.5 2.5 -0.75 3.5 -1 2]; %/
     axis(va); %/
     disp('>> % 	(a more distant view)'');');
     % 	(a more distant view)');
     GAprompt; %/
     disp('>> DrawBivector(v,1.5*v+w); ');
     DrawBivector(v,1.5*v+w); 
     GAtext(0.5*v+0.7*(1.5*v+w)-0.2*unit((v^w)*I3),'v \wedge w');  %/
     axis(va); %/
     GAprompt; %/
     disp('>> DrawBivector(v,-1.5*v+w); ');
     DrawBivector(v,-1.5*v+w); 
     GAtext(0.5*v+0.7*(-1.5*v+w)-0.2*unit((v^w)*I3),'v \wedge w'); %/
     axis(va); %/
     GAprompt; %/
     GAorbiter(360,10); GAprompt; %/
     disp('>> %	(Note that all points spanning same bivector with v');
     %	(Note that all points spanning same bivector with v
     disp('>> %	are on a line parallel to v)');
     %	are on a line parallel to v)
     disp('>> %');
     %
     GAprompt; %/
     disp('>> % 	If no spanning vectors are known ');
     % 	If no spanning vectors are known 
     disp('>> %	we draw the bivector as a disk:');
     %	we draw the bivector as a disk:
     fprintf(1,'\n');     disp('>> draw(v^w,''g''); ');
     draw(v^w,'g'); 
     GAtext(-0.2*unit((v^w)*I3)-w/4,'v \wedge w'); %/ 
     axis(va); %/
     GAprompt; %/
     GAorbiter(360,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 3 ) 
     GAps = '>>>> ';
     disp('>> % 	SPANNING A TRIVECTOR');
     % 	SPANNING A TRIVECTOR
     GAfigure; clc; %/
     disp('>> % 	SPANNING A TRIVECTOR');
     % 	SPANNING A TRIVECTOR
     disp('>> %');
     %
     global dummy; %/
     clf; %/
     disp('>> % 	We can also span a TRIVECTOR:');
     % 	We can also span a TRIVECTOR:
     fprintf(1,'\n');     fprintf(1,'>> e1^e2^e3  ');
     input('');
     e1^e2^e3 %w
     disp('>> % I3  is called the PSEUDOSCALAR of 3-space ');
     % I3  is called the PSEUDOSCALAR of 3-space 
     GAprompt; %/
     disp('>> % 	Again, there is anti-symmetry:');
     % 	Again, there is anti-symmetry:
     fprintf(1,'\n');     fprintf(1,'>> e2^e1^e3  ');
     input('');
     e2^e1^e3 %w
     GAprompt; %/
     disp('>> % 	Here''s how we draw it:');
     % 	Here's how we draw it:
     disp('>> draw(I3);');
     draw(I3);
     GAtext(0.7*e1+0.2*e3,'I3'); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 4 ) 
     GAps = '>>>> ';
     disp('>> % 	SPANNING AND CONTAINMENT');
     % 	SPANNING AND CONTAINMENT
     GAfigure; clc; %/
     disp('>> % 	SPANNING AND CONTAINMENT');
     % 	SPANNING AND CONTAINMENT
     disp('>> %');
     %
     global x B u; %/
     clf; %/
     v = e1; %/
     w = e2; %/
     B = unit(v^w); %/
     draw(B, 'y'); %/
     GAtext(-0.5*v+0.1*B/I3,'B'); %/
     u = B/I3; %/
     draw(u,'r'); %/
     ut = GAtext(0.7*u-0.08*v,'u'); %/
     axis([-0.5 1 -0.5 0.5 -1 1]); %/
     axis off; %/
     disp('>> % We span the trivector u^B, ');
     % We span the trivector u^B, 
     disp('>> % and observe its sign and magnitude as u rotates.');
     % and observe its sign and magnitude as u rotates.
     disp('>> % (click on figure)');
     % (click on figure)
     GAmouse = 1; GAprompt; %/ 
     for i = 0:pi/8:pi %/
     	x = u*cos(i)+v*sin(i); %/
        T = x^B; %/
        title(['u \wedge B = ' num2str(T/I3) '*I_3']);	%/
        delete(ut); %/
        ut = GAtext(0.7*x+0.08*x/(u^v)*norm(u^v),'u'); %/
        if (T/I3 > 1e-16) %/
           draw(x,'r'); %/
       	else if (T/I3 < -1e-16) %/
              draw(x,'b'); %/
        else %/
     		draw(x,'g'); %/
     	end %/
     	end %/
     	axis([-0.5 1 -0.5 0.5 -1 1]); GAprompt; %/ 
     end %/
     GAmouse = 0; %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 5 ) 
     GAps = '>>>> ';
     disp('>> % 	INNER PRODUCT');
     % 	INNER PRODUCT
     GAfigure; clc; %/
     disp('>> % 	INNER PRODUCT');
     % 	INNER PRODUCT
     disp('>> %');
     %
     global x B v w; %/
     clf; %/
     disp('>> % 	The inner product of two vectors is a scalar ');
     % 	The inner product of two vectors is a scalar 
     disp('>> % 	(its magnitude denotes the relative angle)');
     % 	(its magnitude denotes the relative angle)
     fprintf(1,'\n');     disp('>> v = e1;');
     v = e1;
     disp('>> w = unit(e1 + e2);');
     w = unit(e1 + e2);
     draw(v,'b'); %/
     GAtext(0.7*v - 0.1*unit(e2^v)/I3,'v'); %/
     draw(w,'g'); %/
     GAtext(0.7*w + 0.1*unit(e2^w)/I3,'w'); %/
     fprintf(1,'>> inner(v,w)  ');
     input('');
     inner(v,w) %w
     disp('>> % Geometrically, this is a zero-dimensional subspace:');
     % Geometrically, this is a zero-dimensional subspace:
     disp('>> % a weighted point at the origin.');
     % a weighted point at the origin.
     draw(inner(v,w),'r'); %/
     draw(I3/100000,'r'); %/
     GAprompt; %/
     fprintf(1,'\n');     disp('>> % 	The inner product x.B of a vector x and a bivector B is ');
     % 	The inner product x.B of a vector x and a bivector B is 
     disp('>> % 	a vector in B, perpendicular to x');
     % 	a vector in B, perpendicular to x
     disp('>> % 	(its magnitude denotes the relative angle)');
     % 	(its magnitude denotes the relative angle)
     fprintf(1,'\n');     clf; %/
     disp('>> x = unit(e1 + e3); ');
     x = unit(e1 + e3); 
     draw(x,'b'); %/
     GAtext(0.7*x + 0.1*unit(e2^x)/I3,'x'); %/
     disp('>> B = e1^e2; ');
     B = e1^e2; 
     draw(B,'y'); %/
     GAtext(-0.5*v+0.1*B/I3,'B'); %/
     fprintf(1,'>> inner(x,B)  ');
     input('');
     inner(x,B) %w
     axis off; %/
     draw(inner(x,B),'r'); %/
     GAtext(inner(x,B)+0.1*B/I3,'x \bullet B'); GAprompt; %/ 
     disp('>> %	Spinning shows that x.B is "the part of B least like x".');
     %	Spinning shows that x.B is "the part of B least like x".
     GAorbiter(-360,10); %/
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 6 ) 
     GAps = '>>>> ';
     disp('>> % 	DUAL: INNER PRODUCT WITH I3 ');
     % 	DUAL: INNER PRODUCT WITH I3 
     GAfigure; clc; %/
     disp('>> % 	DUAL: INNER PRODUCT WITH I3 ');
     % 	DUAL: INNER PRODUCT WITH I3 
     fprintf(1,'\n');     global x B b; %/
     clf; %/
     disp('>> % 	The inner product with I3 gives the DUAL');
     % 	The inner product with I3 gives the DUAL
     disp('>> %	(the part of I3 least like x).');
     %	(the part of I3 least like x).
     fprintf(1,'\n');     disp('>> x = e1/2 - e2 +e3/3');
     x = e1/2 - e2 +e3/3
     random = unit(pi*e1 + pi/exp(0)*e2 + exp(0)*e3); %/
     fprintf(1,'>> B = inner(x,I3)  ');
     input('');
     B = inner(x,I3) %w
     draw(x,'b'); %/
     GAtext(0.7*x + 0.1*unit(random^x)/I3,'x','b'); %/
     axis off; %/
     axis([-1 1 -1 1 -1 1]); %/
     GAprompt; %/
     disp('>> % 	(The bivector coefficients match the vector coefficients.)');
     % 	(The bivector coefficients match the vector coefficients.)
     fprintf(1,'\n');     draw(B,'g'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAtext(0.3*unit(grade(inner(random,B),1))+0.1*B/I3,'x \bullet I_3'); %/
     GAprompt; %/ 
     GAorbiter(360,10); GAprompt; %/ 
     disp('>> % 	Taking the dual of the bivector gives a vector:');
     % 	Taking the dual of the bivector gives a vector:
     fprintf(1,'\n');     fprintf(1,'>> b = inner(B,I3)  ');
     input('');
     b = inner(B,I3) %w
     draw(b,'r'); %/
     axis([-1 1 -1 1 -1 1]); %/
     GAtext(0.9*b - 0.15*unit(random^b)/I3,'(x \bullet I_3) \bullet I_3','r'); %/ 
     disp('>> ');
     
     disp(' ');     disp('End of GAD sequence.  Returning to Matlab.');
elseif ( GAn == 7 ) 
     GAps = '>>>> ';
     disp('>> % 	GEOMETRIC PRODUCT');
     % 	GEOMETRIC PRODUCT
     GAfigure; clc; %/
     disp('>> % 	GEOMETRIC PRODUCT');
     % 	GEOMETRIC PRODUCT
     disp('>> %');
     %
     global a b; %/
     clf; %/
     disp('>> a = e1 + e3;');
     a = e1 + e3;
     disp('>> b = e1 + e2; ');
     b = e1 + e2; %%
     GAprompt;
     disp('>> %	We use * to denote the geometric product in GABLE:');
     %	We use * to denote the geometric product in GABLE:
     fprintf(1,'>> a*b   ');
     input('');
     a*b %w 
     disp('>> %	Note: scalar + bivector ! ');
     %	Note: scalar + bivector ! 
     GAprompt; %/
     disp('>> %	Order matters!');
     %	Order matters!
     fprintf(1,'>> b*a     ');
     input('');
     b*a    %w
     GAprompt; %/
     disp('>> %	Square of a vector is scalar');
     %	Square of a vector is scalar
     fprintf(1,'>> b*b     ');
     input('');
     b*b    %w
     disp('>> %	Square of a bivector is negative');
     %	Square of a bivector is negative
     fprintf(1,'>> (e1^e2)*(e1^e2)     ');
     input('');
     (e1^e2)*(e1^e2)    %w
     disp('>> %	Every non-null vector has an inverse');
     %	Every non-null vector has an inverse
     fprintf(1,'>> 1/b     ');
     input('');
     1/b    %w
     fprintf(1,'>> b*(1/b)     ');
     input('');
     b*(1/b)    %w
     GAprompt; %/
     disp('>> % 	Inverse formula: ');
     % 	Inverse formula: 
     fprintf(1,'>> b/(b*b)     ');