习题 1. 设 p1=[a1;b1;c1;d1] 和 p2=[a2;b2;c2;d2] 分别表示两个四元数 p1=a1+b1*i+c1*j+d1*k 和 p2=a2+b2*i+c2*j+d2*k ，编写程序求 p=q1*q2。
提示：（1）设 p=a+b*i+c*j+d*k 则可用 p=[a; b; c; d] 表示.
   （2）令
          A=[a1 -b1 -c1 -d1
            b1 a1 -d1 c1
            c1 d1 a1 -b1
            d1 -c1 b1 a1],
则 p=A*p2.

分析与解答：

这是求解四元数乘法的另一种方式，把p1扩充成矩阵的形式，p2写成向量的形式，p=p1*p2就是一般意义下的矩阵乘法。p1的每一行与p2的 列向量的乘积分别对应p的常数项,i, j, k.所以矩阵p1的每一行是按着四元数项城的规则排列的。输入p1,调用函数 A2111，使p1按着常数项,i, j, k.的顺序排列。

代码;

function  p=timea(B,C)

[B0,B1,B2,B3]=A2111(B);

% [C0,C1,C2,C3]=A2A0123(C);

A0=[B0,-B1,-B2,-B3];

A1=[B1,B0,-B3,B2];

A2=[B2,B3,B0,-B1];

A3=[B3,-B2,B1,B0];

A=[A0 ;A1; A2 ;A3];

p=A*C;

function [A0,A1,A2,A3]=A2111(A)

n=size(A,2)/4;

A0=A(:,1:n);

A1=A(:,(n+1):(2*n));

A2=A(:,(2*n+1):(3*n));

A3=A(:,(3*n+1):(4*n))

试验结果：

>> clear

>> syms a1 b1 c1 d1 a2 b2 c2 d2

>> B=[a1 b1 c1 d1];

>> C=[a2; b2; c2 ;d2];

>> timea(B,C)

ans =

 a1*a2-b1*b2-c1*c2-d1*d2

 b1*a2+a1*b2-d1*c2+c1*d2

 c1*a2+d1*b2+a1*c2-b1*d2

 d1*a2-c1*b2+b1*c2+a1*d2


习题 2. 设 q1= 2+1*i+3*j+4*k 和 q2=0+1*i+1*j-1*k，利用习题1中编写的程序计算 p=p1*p2, 并且利用 Bj Bi Bk 画图解释四元数乘积与旋转之间的关系。

分析与解答：

对于p=p1*p2中，如果p2的常数项为0，p1的后三行元素I,j,k对应的I,j,k形成矩阵B.并且Bj Bi Bk求和形成矩阵B，其中 Bj=[a1/2,0,c1;0,0,0;-c1,0,a1/2]; Bk=[a1/2,-d1,0;d1,a1/2,0;0,0,0]; Bi=[0,0,0;0,a1/2,-b1;0,b1,a1/2];

Bj作用于p2列向量，即绕着Bj旋转， Bi Bk类似，形成向量cj,ck,ci,经过一系列的旋转变化p=Bj*cj+ Bk*ck+ Bi*ci。

代码：

A = [1,1,0;

    -3 ,5, 0 ];    %输入变化前后的向量c

[mm,nn] = size(A);

hold on

grid on

    plot3(A(:,1),A(:,2),A(:,3),'o')  

for i = 1:mm

    X1 = [A(i,1);0];

    Y1 = [A(i,2);0];

    Z1 = [A(i,3);0];

    plot3(X1,Y1,Z1)

end

view(-20,40)

xlabel( 'x');

ylabel( 'y');

zlabel( 'z');

试验结果：

调用函数timea(B,C)，计算p=p1*p2，画出Bj Bi Bk作用下p的形成过程。

A=[2 -1 -3 -4; 1  2  -4  3;3  4  2  -1;4 -3 1 2 ];

>> clear

>> B=[2 1 3 4];

>> C=[0;1;1;-1];

>> timea(B,C)

ans =

     0

    -5

     7

-4

Bj=[1,0,3;0,0,0;-3,0,1];  cj=[1,0,-1]'; c=[-2 ,0 ,-4]

两条线分别表示在Bj作用前后变化的向量。 

Bk=[1,-4,0;4,1,0;0,0,0];  ck=[1,1,0]'; c=[ -3 ,5, 0]

Bi=[0,0,0;0,1,-1;0,1,1];  ci=[0,1,-1]'; c=[0, 2 ,0]