Calculator


合并了《Calculator V1.8.6.0（2013-10-31 10:33）》、《山寨版Matlab——Calculator V5.5beta(5.5版本更新 By2013/11/21 )2013-12-4 22:37》这几篇QQ日志。

如果A是一个m*n的矩阵，而B是一个p*q的矩阵，克罗内克积AⓧB则是一个mp*nq的分块矩阵。

矩阵的克罗内克积对应于线性映射的抽象张量积。并矢积是克罗内克积的特殊情况。克罗内克积是张量积的特殊形式，因此满足双线性与结合律，不符合交换律。

3行2列矩阵a与2行2列矩阵b的克罗内克积为6行4列的矩阵c

2行2列矩阵a与2行2列矩阵的克罗内克积为4行4列的矩阵c

如果A（m行n列），B（p行q列），C（n行c列），D（q行d列）是4个矩阵，且矩阵乘积AC（m行c列）和BD（p行d列）存在，那么：(AⓧB) (CⓧD)=
ACⓧBD=mp行cd列的矩阵=mp行nq列矩阵乘以nq行cd列矩阵

混合积的性质：它混合了通常的矩阵乘积和克罗内克积

AⓧB是可逆的当且仅当A（m行m列）和B（p行p列）是可逆的，其逆矩阵为：(AⓧB)^(-1)=A^(-1)ⓧB^(-1) =（mp行mp列）

如果A是n行n列矩阵，B是m行m列矩阵，I_k表示k行k列单位矩阵，那么我们可以定义克罗内克和⊕为：A⊕B=AⓧI_m+I_nⓧB=nm行nm列的矩阵

克罗内克积转置运算符合分配律：

如果A（m行n列），B（p行q列），(AⓧB)^(T)=A^(T)ⓧB^(T)=(mp行nq列)^(T)=nq行mp列的矩阵

rank(AⓧB)=rank(A)+rank(B)
功能:对矩阵进行Kronecker张量积.

格式:Kron(a,b)

说明:如果a是m*n的,b是p*q的,则kron(a,b)是大小(m*p)*(n*q)的矩阵

例子:

a=[1 2

3 1]

a =

[ 1.00000000000000 2.00000000000000

3.00000000000000 1.00000000000000 ]

b=[0 3

2 1]

b =

[ 0.00000000000000 3.00000000000000

2.00000000000000 1.00000000000000 ]

c=Kron(a,b)

c =

[ 0.00000000000000 3.00000000000000 0.00000000000000 6.00000000000000

2.00000000000000 1.00000000000000 4.00000000000000 2.00000000000000

0.00000000000000 9.00000000000000 0.00000000000000 3.00000000000000

6.00000000000000 3.00000000000000 2.00000000000000 1.00000000000000 ]

a=[1

2

3]

a =

[ 1.00000000000000

2.00000000000000

3.00000000000000 ]

b=[4 5 6]

b =

[ 4.00000000000000 5.00000000000000 6.00000000000000 ]

c=Kron(a,b)

c =

[ 4.00000000000000 5.00000000000000 6.00000000000000

8.00000000000000 10.0000000000000 12.0000000000000

12.0000000000000 15.0000000000000 18.0000000000000 ]

d=mul(a,b)

d =

[ 4.00000000000000 5.00000000000000 6.00000000000000

8.00000000000000 10.0000000000000 12.0000000000000

12.0000000000000 15.0000000000000 18.0000000000000 ]

功能:求矩阵的2范数.谱范数.

格式:Cond2(a)

说明:a是任意维数的矩阵.本函数执行原理是,返回对称矩阵a'*a的最大特征值的平方根.其中a'是a的转置.

a=[1 -2

-3 4]

a =

[ 1.00000000000000 -2.00000000000000

-3.00000000000000 4.00000000000000 ]

cond2(a)

ans =

[ 5.46498570421904 ]

矩阵A的2范数就是 A的转置乘以A矩阵特征根 最大值的开根号

如A={ 1 -2

-3 4 }

那么A的2范数就是msgbox (15+221^0.5)^0.5 '=5.46498570421904

范数是一种实泛函，满足非负性、齐次性、三角不等式

要理解矩阵的算子范数，首先要理解向量范数的内涵。矩阵的算子范数，是由向量范数导出的。

范数理论是矩阵分析的基础，度量向量之间的距离、求极限等都会用到范数，范数还在机器学习、模式识别领域有着广泛的应用。

矩阵范数反映了线性映射把一个向量映射为另一个向量，向量的"长度"缩放的比例。

由矩阵算子范数的定义形式可知，矩阵A把向量x映射成向量Ax，取其在向量x范数为1所构成的闭集下的向量Ax范数最大值作为矩阵A的范数，即矩阵对向量缩放的比例的上界，矩阵的算子范数是相容的。由几何意义可知，矩阵的算子范数必然大于等于矩阵谱半径（最大特征值的绝对值），矩阵算子范数对应一个取到向量Ax范数最大时的向量x方向，谱半径对应最大特征值下的特征向量的方向。而矩阵的奇异值分解SVD，分解成左右各一个酉阵，和拟对角矩阵，可以理解为对向量先作旋转、再缩放、最后再旋转，奇异值，就是缩放的比例，最大奇异值就是谱半径的推广，所以，矩阵算子范数大于等于矩阵的最大奇异值，酉阵在此算子范数的意义下，范数大于等于1。此外，不同的矩阵范数是等价的。

功能:方阵求逆,其采用的是LU分解算法

格式:inv(a)//a为方阵变量

a=[0.2368 0.2471 0.2568 1.2671

1.1161 0.1254 0.1397 0.149

0.1582 1.1675 0.1768 0.1871

0.1968 0.2071 1.2168 0.2271]

a =

[ 0.23680000000000 0.24710000000000 0.25680000000000 1.26710000000000

1.11610000000000 0.12540000000000 0.13970000000000 0.14900000000000

0.15820000000000 1.16750000000000 0.17680000000000 0.18710000000000

0.19680000000000 0.20710000000000 1.21680000000000 0.22710000000000 ]

b=inv(a)

b =

[ -0.0859207504780 0.93794426823404 -0.0684372042645 -0.0796077151837

-0.1055899132073 -0.0885243235004 0.90598255638825 -0.0991908105397

-0.1270733117900 -0.1113511370480 -0.1169667064884 0.87842529094384

0.85160581464323 -0.1354556628418 -0.1401825503018 -0.1438074804470 ]

ublas求矩阵的逆的输出及源代码：

mtxA矩阵[4,4]((0.2368,0.2471,0.2568,1.2671),(1.1161,0.1254,0.1397,0.149),(0.1582

,1.1675,0.1768,0.1871),(0.1968,0.2071,1.2168,0.2271))

mtxA逆矩阵[4,4]((-0.0859208,0.937944,-0.0684372,-0.0796077),(-0.10559,-0.0885243

,0.905983,-0.0991908),(-0.127073,-0.111351,-0.116967,0.878425),(0.851606,-0.1354

56,-0.140183,-0.143807))

mtxA逆矩阵的逆矩阵 [4,4]((0.2368,0.2471,0.2568,1.2671),(1.1161,0.1254,0.1397,0.1

49),(0.1582,1.1675,0.1768,0.1871),(0.1968,0.2071,1.2168,0.2271))

// REMEMBER to update "lu.hpp" header includes from boost-CVS

#include <boost/numeric/ublas/vector.hpp>

#include <boost/numeric/ublas/vector_proxy.hpp>

#include <boost/numeric/ublas/matrix.hpp>

#include <boost/numeric/ublas/triangular.hpp>

#include <boost/numeric/ublas/lu.hpp>

#include <boost/numeric/ublas/io.hpp>

#include <boost/numeric/ublas/vector.hpp>
#include <iostream>
namespace ublas = boost::numeric::ublas;
//! \brief 求矩阵的逆

//! \param input[in] 被求矩阵

//! \param inverse[out] 逆矩阵

/* Matrix inversion routine.

Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */

template<class T>

bool InvertMatrix (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {

using namespace boost::numeric::ublas;

typedef permutation_matrix<std::size_t> pmatrix;

// create a working copy of the input

matrix<T> A(input);

// create a permutation matrix for the LU-factorization

pmatrix pm(A.size1());
// perform LU-factorization

int res = lu_factorize(A,pm);

if( res != 0 ) return false;
// create identity matrix of "inverse"

inverse.assign(ublas::identity_matrix<T>(A.size1()));
// backsubstitute to get the inverse

lu_substitute(A, pm, inverse);
return true;

}

int main()

{

namespace ublas = boost::numeric::ublas;

ublas::matrix<double> mTmp(4,4);

/*

' 原矩阵

mtxA(1, 1) = 0.2368: mtxA(1, 2) = 0.2471: mtxA(1, 3) = 0.2568: mtxA(1, 4) = 1.2671

mtxA(2, 1) = 1.1161: mtxA(2, 2) = 0.1254: mtxA(2, 3) = 0.1397: mtxA(2, 4) = 0.149

mtxA(3, 1) = 0.1582: mtxA(3, 2) = 1.1675: mtxA(3, 3) = 0.1768: mtxA(3, 4) = 0.1871

mtxA(4, 1) = 0.1968: mtxA(4, 2) = 0.2071: mtxA(4, 3) = 1.2168: mtxA(4, 4) = 0.2271

*/
#if 1

//! \brief 求逆矩阵

ublas::matrix<double> mtxA(4,4), MAInv(4,4);

mtxA(0, 0) = 0.2368; mtxA(0, 1) = 0.2471; mtxA(0, 2) = 0.2568; mtxA(0, 3) = 1.2671;

mtxA(1, 0) = 1.1161; mtxA(1, 1) = 0.1254; mtxA(1, 2) = 0.1397; mtxA(1, 3) = 0.149;

mtxA(2, 0) = 0.1582; mtxA(2, 1) = 1.1675; mtxA(2, 2) = 0.1768; mtxA(2, 3) = 0.1871;

mtxA(3, 0) = 0.1968; mtxA(3, 1) = 0.2071; mtxA(3, 2) = 1.2168; mtxA(3, 3) = 0.2271;

std::cout <<"mtxA矩阵"<< mtxA <<std::endl;

InvertMatrix(mtxA, MAInv);

InvertMatrix(mtxA, mTmp);

std::cout <<"mtxA逆矩阵"<< mTmp << std::endl;

InvertMatrix(mTmp, mTmp);

std::cout <<"mtxA逆矩阵的逆矩阵 " <<mTmp <<std::endl;

std::cout <<"" <<std::endl;

#endif
system("pause");

return 0;

}

功能:矩阵相乘

格式:mul(a,b)

a=[1 3 -2 0 4

-2 -1 5 -7 2

0 8 4 1 -5

3 -3 2 -4 1]

a =

[ 1.00000000000000 3.00000000000000 -2.0000000000000 0.00000000000000 4.00000000000000

-2.0000000000000 -1.0000000000000 5.00000000000000 -7.0000000000000 2.00000000000000

0.00000000000000 8.00000000000000 4.00000000000000 1.00000000000000 -5.0000000000000

3.00000000000000 -3.0000000000000 2.00000000000000 -4.0000000000000 1.00000000000000 ]

b=[4 5 -1

2 -2 6

7 8 1

0 3 -5

9 8 -6]

b =

[ 4.00000000000000 5.00000000000000 -1.0000000000000

2.00000000000000 -2.0000000000000 6.00000000000000

7.00000000000000 8.00000000000000 1.00000000000000

0.00000000000000 3.00000000000000 -5.0000000000000

9.00000000000000 8.00000000000000 -6.0000000000000 ]

c=mul(a,b)

c =

[ 32.0000000000000 15.0000000000000 -9.0000000000000

43.0000000000000 27.0000000000000 24.0000000000000

-1.0000000000000 -21.000000000000 77.0000000000000

29.0000000000000 33.0000000000000 -5.0000000000000 ]

相应的C++输出及源代码（基于boost::multi_array的矩阵相乘）：

#include <iostream>

#include "boost/multi_array.hpp"

using namespace std;
typedef boost::multi_array<int, 2> matrix;
matrix matrix_multiply(matrix& a,matrix& b)

{

matrix::index row=a.shape()[0];

matrix::index col=b.shape()[1];

matrix c(boost::extents[row][col]);
for (matrix::index i=0; i!=a.shape()[0]; ++i)

for (matrix::index j=0; j!=b.shape()[1]; ++j)

for (matrix::index k=0; k!=a.shape()[1]; ++k)

c[i][j]+=a[i][k]*b[k][j];
return c;

}
void print(const matrix& m)

{

for (matrix::index i=0; i!=m.shape()[0]; cout<<endl,++i)

for (matrix::index j=0; j!=m.shape()[1]; ++j)

cout<<m[i][j]<<" ";

}
int main() {
int valuesA[] = {

1, 3, -2,0,4,

-2,-1, 5,-7,2,

0,8,4,1,-5,

3,-3,2,-4,1

};

int valuesB[] = {

4, 5,-1,

2, -2,6,

7,8,1,

0,3,-5,

9,8,-6

};

const int valuesA_size = 20;

const int valuesB_size = 15;

matrix A(boost::extents[4][5]);

matrix B(boost::extents[5][3]);

A.assign(valuesA,valuesA + valuesA_size);

B.assign(valuesB,valuesB + valuesB_size);
/*

Dim mtxA(4, 5) As Double

Dim mtxB(5, 3) As Double

Dim mtxC(4, 3) As Double

mtxA(1, 1) = 1: mtxA(1, 2) = 3: mtxA(1, 3) = -2: mtxA(1, 4) = 0: mtxA(1, 5) = 4

mtxA(2, 1) = -2: mtxA(2, 2) = -1: mtxA(2, 3) = 5: mtxA(2, 4) = -7: mtxA(2, 5) = 2

mtxA(3, 1) = 0: mtxA(3, 2) = 8: mtxA(3, 3) = 4: mtxA(3, 4) = 1: mtxA(3, 5) = -5

mtxA(4, 1) = 3: mtxA(4, 2) = -3: mtxA(4, 3) = 2: mtxA(4, 4) = -4: mtxA(4, 5) = 1
mtxB(1, 1) = 4: mtxB(1, 2) = 5: mtxB(1, 3) = -1

mtxB(2, 1) = 2: mtxB(2, 2) = -2: mtxB(2, 3) = 6

mtxB(3, 1) = 7: mtxB(3, 2) = 8: mtxB(3, 3) = 1

mtxB(4, 1) = 0: mtxB(4, 2) = 3: mtxB(4, 3) = -5

mtxB(5, 1) = 9: mtxB(5, 2) = 8: mtxB(5, 3) = -6

*/
cout<<"matrix A"<<endl;

print(A);

cout<<endl;cout<<"*"<<endl;cout<<"matrix B"<<endl;

print(B);

cout<<endl;cout<<"="<<endl;cout<<"matrix C"<<endl;

print(matrix_multiply(A,B));

cout<<endl;
system("pause");

return 0;

}
tan(0.25+0.25i)

ans = 0.239090115629191 + 0.259870880485541i

(9+11j)(56+3j)=471+643j

(9+11i)*(56+3i)

ans = 471 + 643i

x=[1 0

0 1]

x =

[ 1.00000000000000 0.00000000000000

0.00000000000000 1.00000000000000 ]

sin(x)

ans =

[ 0.84147098480789 0.00000000000000

0.00000000000000 0.84147098480789 ]
x=[0.25 -0.25

0.25 0.25]

x =

[ 0.25000000000000 -0.2500000000000

0.25000000000000 0.25000000000000 ]

tan(x)

ans =

[ 0.25534192122103 -0.2553419212210

0.25534192122103 0.25534192122103 ]
x=[0.25 0.25

-0.25 0.25]

x =

[ 0.25000000000000 0.25000000000000

-0.2500000000000 0.25000000000000 ]

tan(x)

ans =

[ 0.25534192122103 0.25534192122103

-0.2553419212210 0.25534192122103 ]

CTan(0.25+0.25j)=0.2391+(-0.26j)

功能:采用高斯消元法求解ax=b这类方程.

格式:Gauss(a,b)

注意:a是方阵,b是n*1的矩阵
a=[0.2368 0.2471 0.2568 1.2671

0.1968 0.2071 1.2168 0.2271

0.1581 1.1675 0.1768 0.1871

1.1161 0.1254 0.1397 0.1490]

a =

[ 0.23680000000000 0.24710000000000 0.25680000000000 1.26710000000000

0.19680000000000 0.20710000000000 1.21680000000000 0.22710000000000

0.15810000000000 1.16750000000000 0.17680000000000 0.18710000000000

1.11610000000000 0.12540000000000 0.13970000000000 0.14900000000000 ]

b=[1.8471

1.7471

1.6471

1.5471]

b =

[ 1.84710000000000

1.74710000000000

1.64710000000000

1.54710000000000 ]

x=gauss(a,b)

x =

[ 1.04057667941935

0.98705076839213

0.93504033393356

0.88128232948438 ]

相应的C程序输出及源代码：
#include <stdio.h>

#include <math.h>

#include <windows.h>

extern "C" _declspec(dllexport)int __stdcall agaus(double *a,double *b,int n);

int __stdcall agaus(double *a,double *b,int n)

{

int *js,l,k,i,j,is,p,q;

double d,t;

js=(int *)malloc(n*sizeof(int));

l=1;

for (k=0;k<=n-2;k++)

{ d=0.0;

for (i=k;i<=n-1;i++)

for (j=k;j<=n-1;j++)

{ t=fabs(a[i*n+j]);

if (t>d) { d=t; js[k]=j; is=i;}

}

if (d+1.0==1.0) l=0;

else

{ if (js[k]!=k)

for (i=0;i<=n-1;i++)

{ p=i*n+k; q=i*n+js[k];

t=a[p]; a[p]=a[q]; a[q]=t;

}

if (is!=k)

{ for (j=k;j<=n-1;j++)

{ p=k*n+j; q=is*n+j;

t=a[p]; a[p]=a[q]; a[q]=t;

}

t=b[k]; b[k]=b[is]; b[is]=t;

}

}

if (l==0)

{ free(js); printf("fail\n");

return(0);

}

d=a[k*n+k];

for (j=k+1;j<=n-1;j++)

{ p=k*n+j; a[p]=a[p]/d;}

b[k]=b[k]/d;

for (i=k+1;i<=n-1;i++)

{ for (j=k+1;j<=n-1;j++)

{ p=i*n+j;

a[p]=a[p]-a[i*n+k]*a[k*n+j];

}

b[i]=b[i]-a[i*n+k]*b[k];

}

}

d=a[(n-1)*n+n-1];

if (fabs(d)+1.0==1.0)

{ free(js); printf("fail\n");

return(0);

}

b[n-1]=b[n-1]/d;

for (i=n-2;i>=0;i--)

{ t=0.0;

for (j=i+1;j<=n-1;j++)

t=t+a[i*n+j]*b[j];

b[i]=b[i]-t;

}

js[n-1]=n-1;

for (k=n-1;k>=0;k--)

if (js[k]!=k)

{ t=b[k]; b[k]=b[js[k]]; b[js[k]]=t;}

free(js);

return(1);

}
int main()

{

typedef int (__stdcall *Func)(double *a,double *b,int n);

//HINSTANCE hDLL;

int i;

Func f=agaus;

/*

x(0)=1.040577e+000

x(1)=9.870508e-001

x(2)=9.350403e-001

x(3)=8.812823e-001

*/

static double a[4][4]={ {0.2368,0.2471,0.2568,1.2671},

{0.1968,0.2071,1.2168,0.2271},

{0.1581,1.1675,0.1768,0.1871},

{1.1161,0.1254,0.1397,0.1490} };

static double b[4]={1.8471,1.7471,1.6471,1.5471};

//hDLL=LoadLibrary("mathlib.dll");

//f=(Func)GetProcAddress(hDLL, "_agaus@12");

if (f(&a[0][0],&b[0],4)!=0)

{

for (i=0;i<=3;i++)

printf("x(%d)=%e\n",i,b[i]);

}

//FreeLibrary(hDLL);

getchar();

return 0;

}

数据处理之傅里叶变换Dft和傅里叶逆变换IDft：

a=[1

2

3

4]

a =

[ 1.00000000000000

2.00000000000000

3.00000000000000

4.00000000000000 ]

b=dft(a,4)

b =

[ 10.0000000000000 0.00000000000000

-2.0000000000000 2.00000000000000

-2.0000000000000 -9.796850830E-16

-2.0000000000000 -2.0000000000000 ]

c=idft(b)

c =

[ 1.00000000000000 -5.551115123E-16

2.00000000000000 -2.775557561E-16

3.00000000000000 -1.110223024E-16

4.00000000000000 2.7755575615E-16 ]

相应的C程序输出及源代码：

请输入N：4

1

2

3

4

test and verify:

10.000000

-2.000000

-2.000000

-2.000000
0.000000

2.000000

-0.000000

-2.000000

1.000000

2.000000

3.000000

4.000000
-0.000000

0.000000

-0.000000

-0.000000

#include<windows.h>

#include<stdio.h>

#include<stdlib.h>

#include<math.h>
//1D-DFT.exe利用定义计算任意点实序列的一维DFT

typedef int(*pDFT1)(double pr[],double pi[],double fr[],double fi[],int N,int inv);
int DFT1(double pr[],double pi[],double fr[],double fi[],int N,int inv)

{

double sr,si;

if(inv==1)

{

for(int j=0;j<N;j++)

{

sr=0.0;

si=0.0;

for(int k=0;k<N;k++)

{

sr+=pr[k]*cos(-6.2831853068*j*k/N);

si+=pr[k]*sin(-6.2831853068*j*k/N);

}

fr[j]=sr;

fi[j]=si;

}
for(int i=0;i<N; i++)

{

printf("%lf\n",fr[i]);

}

printf("\n");

for(int i=0; i<N; i++)

{

printf("%lf\n",fi[i]);

}

printf("\n\n");
}
if(inv==-1)

{

for(int k=0;k<N;k++)

{

sr=0.0;

si=0.0;

for(int j=0;j<N;j++)

{

sr+= fr[j]*cos(6.2831853068*j*k/N)-fi[j]*sin(6.2831853068*j*k/N);

si+= fr[k]*sin(6.2831853068*j*k/N)+fi[j]*cos(6.2831853068*j*k/N);

}

pr[k]=sr/N;

pi[k]=si/N;

}
for(int i=0;i<N; i++)

{

printf("%lf\n",pr[i]);

}

printf("\n");

for(int i=0; i<N; i++)

{

printf("%lf\n",pi[i]);

}

printf("\n\n");

}

return 0;

}
int main(void)

{

#if 1

int i,j,k,N;

double *pr,*pi,*fr,*fi,sr,si;

printf("请输入N：");

scanf("%d",&N);

pr=(double *)malloc(N*sizeof(double));

pi=(double *)malloc(N*sizeof(double));

fr=(double *)malloc(N*sizeof(double));

fi=(double *)malloc(N*sizeof(double));

for (int i=0; i<N;i++)

{

scanf("%lf",pr+i);

*(pi+i)=0;

}

printf("test and verify:\n");

DFT1(pr,pi,fr,fi,N,1);

DFT1(pr,pi,fr,fi,N,-1);

free(pr);

pr=NULL;

free(pi);

pi=NULL;

free(fr);

fr=NULL;

free(fi);

fi=NULL;

#endif
#if 0

//_FFT1@16一维光学DFT的计算（调用FFT1）

typedef int(*pFFT1)(float a_rl[], float a_im[], int ex,int inv);

/*

请输入N,K：4

2

1

2

3

4

-1.000000,0.000000

1.000000,1.000000

5.000000,0.000000

1.000000,-1.000000

1.000000,0.000000

2.000000,-0.000000

3.000000,0.000000

4.000000,0.000000

*/
HINSTANCE hDLL=NULL;

pFFT1 FFT1=NULL;

float *p1=NULL,*p2=NULL;

int N,K;

printf("请输入N,K：");

scanf("%d%d",&N,&K);

p1=(float *)malloc(N*sizeof(float));

p2=(float *)malloc(N*sizeof(float));

for(int i=0;i<N;i++)

{

scanf("%f",p1+i);

*(p2+i)=0;

}

hDLL=LoadLibrary("dft.dll");

FFT1=(pFFT1)GetProcAddress(hDLL, "FFT1");

int a =FFT1(p1,p2,K,1);

for(int i=0;i<N;i++)

{

printf("%f,%f\n",p1[i],p2[i]);

}

a=FFT1(p1,p2,K,-1);

for(int i=0;i<N;i++)

{

printf("%f,%f\n",p1[i],p2[i]);

}

FreeLibrary(hDLL);

free(p1);

p1=NULL;

free(p2);

p2=NULL;

#endif
system("pause");

return 0;

}

正弦积分_lsinn@8

sinint(0)

ans =

[ 0.00000000000000 ]

sinint(1)

ans =

[ 0.94608307041448 ]

sinint(1000)

ans =

[ 1.57023311828152 ]

余弦积分

cosint(0)

ans =

[ -80.013262589890 ]

cosint(0.5)

ans =

[ -0.6777840788097 ]

cosint(1)

ans =

[ -0.6625960771110 ]

cosint(1.5)

ans =

[ -1.0296436828302 ]

cosint(1000)

ans =

[ -999.99917370536 ]

另一种余弦积分_lcoss@8

Ci(0.00)= -80.0132626

Ci(0.50)= -0.1777841

Ci(1.00)= 0.3374039

Ci(1.50)= 0.4703563

Ci(1000.00)= 0.0008189

伽马函数值Γ(0.5)

gamma(0.5)

ans =

[ 1.77245385090559 ]

由伽马函数表示的贝塔函数值B(0.5,0.5)=Γ(0.5)*Γ(0.5)/Γ(1.0))

Beta(0.5,0.5)

ans =

[ 3.14159265359006 ]

功能:求方阵行列式.此方法精度高,但运算的方阵的阶数最好不要超过5.
格式:Det(a)
说明:a为方阵
例子:
a =
[ 24.670817016005 97.0357947969045 89.6194347132088
2.72601312153321 14.5667740211667 58.9581970399982
6.19456609999508 60.2519491967987 8.23347480419719 ]
执行命令
b=det(a)
后输出
b =
[ -44785.8743735347 ]



判断矩阵是否是正交矩阵，是否是旋转矩阵。

a=[0.99770898 0.06753442 0.00398691

-0.06752640 0.99771525 -0.00211390

-0.00412057 0.00183984 0.99998982]

a =

[ 0.99770898000000 0.06753442000000 0.00398691000000

-0.06752640000000 0.99771525000000 -0.00211390000000

-0.00412057000000 0.00183984000000 0.99998982000000 ]

b=det(a)

b =

[ 1.00000000483672 ]

b=inv(a)

b =

[ 0.99770897767089 -0.06752639689650 -0.00412056189788

0.06753442272752 0.99771524683865 0.00183983532449

0.00398691773663 -0.00211390518455 0.99998981582996 ]

a'

ans =

[ 0.99770898000000 -0.06752640000000 -0.00412057000000

0.06753442000000 0.99771525000000 0.00183984000000

0.00398691000000 -0.00211390000000 0.99998982000000 ]

c=a'

c =

[ 0.99770898000000 -0.06752640000000 -0.00412057000000

0.06753442000000 0.99771525000000 0.00183984000000

0.00398691000000 -0.00211390000000 0.99998982000000 ]

IsEqu(b,c,3)

2矩阵相等

20131207单位四元数和三阶幺模矩阵的互转公式实际上给出了拓扑学中S^3和SO(3)XI_2同胚的构造性证明。

有限维可交换Lie群：

一维可交换Lie群——用矩阵（线性变换）的形式来表示绕固定轴的旋转群SO(2)：

{{x'},{y'}}={{cosθ,-sinθ},{sinθ,cosθ}}{{x},{y}}={{xcosθ-ysinθ},{xcosθ-ysinθ}}

----这里的θ就是单位复数cosθ+isinθ的辐角主值，即θ=arg(cosθ+isinθ)=atan2(sinθ,cosθ)

----虚数i代表[逆时针]旋转90度，i*i=-1代表[逆时针]旋转180度，(-1)(-1)=1代表逆时针旋转360度，当然转回来了。

【

复数的矩阵表示：

1.令a+bi对应{{a,b},{-b,a}}，记σ(a+bi)={{a,b},{-b,a}}，例如σ(cosα+isinα)={{cosα,sinα},{-sinα,cosα}}/*顺时针旋转辐角主值α*/,α=arg(a+bi)=atan2(b,a)

2.令a+bi对应{{a,-b},{b,a}}，记σ(a+bi)={{a,-b},{b,a}}，例如σ(cosθ+isinθ)={{cosθ,-sinθ},{-sinθ,cosθ}}/*逆时针旋转辐角主值θ*/,θ=arg(a+bi)=atan2(b,a)

】

一个三维紧致非单连通不可交换单Lie群——三维转动群SO(3)也称为特殊正交群.三维空间绕固定点的一个转动，习惯用Euler角（θ，Φ，Ψ）来描述一个转动。不是单连通的流形。

g^Ψ_z表示绕z轴转Ψ角

g^θ_x表示绕x轴转θ角

g^Φ_y表示绕y轴转Φ角

于是g=g^Φ_yg^θ_xg^Ψ_z

令θ，Φ=0，即绕固定点O和固定轴z轴旋转Ψ

则{{x'},{y'}}={{cosΨ,-sinΨ},{sinΨ,cosΨ}}{{x},{y}}={{xcosΨ-ysinΨ},{xcosΨ-ysinΨ}}

{{x'},{y'},{z'}}={{cosΨ,-sinΨ,0},{sinΨ,cosΨ,0},{0,0,1}}{{x},{y},{z}}={{xcosΨ-ysinΨ},{xcosΨ-ysinΨ},{z}}

仅满足g^tg=gg^t[而不要求det g>0]的线性变换所构成的群称为O(3)——正交群。

【

三维紧致单连通单Lie群SU(2)[幺正幺模]的表示：

单位四元数组成Lie群SU(2)，其基本表示为Pauli矩阵。

SU(2)的基本表示（二维表示）是g={{z_1,-~z_2},{z_2,~z_1}}∈SU(2)

其中复数z_1,z_2满足条件：|z_1|^2+|z_2|^2=1

令z_1=t-iz,z_2=y-ix

则约束可表示为t^2+z^2+y^2+x^2=1

SU(2)群流形相当于R^4中的单位球S^3

|i|=|j|=|k|=1

三维紧致非单连通不可交换单Lie群SO(3)，群元表示为行列式值为1的3×3正交矩阵

思考：单Lie群的定义与SO(3)=SU(2)/Z_2=S^3/Z_2=RP(3)的关系？

SU(2)和SO(3)的Lie代数：L(SO(3))=L(SU(2))=A_1[三维单的半单Lie代数，R^3上附加三维叉乘为李括号的Lie代数]

单位复数a+bi转换到矩阵：

Matrix={{a,-b},{b,a}}

1，i对应的矩阵分别为：

{{1,0},{0,1}}

{{0,-1},{1,0}}

{1,i,-1,-i}=C_4

C_4中的4个群元代表E^2中的4个旋转

单位四元数w+xi+yj+zk转换到矩阵：

Matrix={{1-2y^2-2z^2,2xy-2wz,2xz+2wy},{2xy+2wz,1-2x^2-2z^2,2yz-2wx},{2xz-2wy,2yz+2wx,1-2x^2-2y^2}}

1,i,j,k对应的矩阵分别为：

{{1,0,0},{0,1,0},{0,0,1}},

{{1,0,0},{0,-1,0},{0,0,-1}}

{{-1,0,0},{0,1,0},{0,0,-1}}

{{-1,0,0},{0,-1,0},{0,0,1}}

{1,i,-1,-i,j,k,-j,-k}=Q_8

Q_8中的8个群元代表E^3中的8个旋转

】

三维转动群SO(3)是3维连通紧致单线性Lie群，相应的实Lie代数so(3)是3维1秩紧致实单Lie代数。

二维Lorentz群SO(2，1)是3维非紧致单线性Lie群，相应的实Lie代数so(2，1)是3维1秩非紧致实单Lie代数。

非交换群（非阿贝尔群）的例子有各种李群（是拓扑群也是微分流形）——线性群之有限维不可交换李群：

n阶实[一般]/[全]线性群GL(n,R)扩张为n阶复[一般]/[全]线性群GL(n,C)

n阶实特殊线性群SL(n,R)扩张为n阶复特殊线性群SL(n,C)

该软件目前还没有四元数和旋转矩阵互相转换的功能。

1,i,j,k的欧拉角向量(omiga,phi,kappa)=(i_pi,j_pi,k_pi)分别为(0,0,0)，(pi,0,0)，(0,pi,0)，(0,0,pi)。

RotMat r =

0.600266 0.152151 -0.785195

-0.799692 0.130352 -0.58609

0.013178 0.979724 0.19992

Quat q =

-0.56347 0.287301 0.342528 0.694719

RotMat _r =

0.600266 0.152151 -0.785195

-0.799692 0.130352 -0.58609

0.013178 0.979724 0.19992

struct Quat {

ValType _v[4];//x, y, z, w

//phi = 1.32148229302237 ; omiga = 0.626224465189316 ; kappa = -1.4092143985971;

四元数1={{1,0,0},{0,1,0},{0,0,1}}

phi = 0 ; omiga = 0 ; kappa = 0;

RotMat r =

1 -0 -0

0 1 -0

0 0 1

Quat q =

-0 0 -0 1

RotMat _r =

1 -0 0

0 1 -0

0 0 1

四元数j={{-1,0,0},{0,1,0},{0,0,-1}}

phi = atan2((double)0,-1) ; omiga = 0 ; kappa = 0;

RotMat r =

-1 0 -1.22465e-016

0 1 -0

1.22465e-016 -0 -1

Quat q =

0 1 -0 6.12323e-017

RotMat _r =

-1 0 -1.22465e-016

0 1 0

1.22465e-016 -0 -1

四元数i={{1,0,0},{0,-1,0},{0,0,-1}}

phi = 0 ; omiga = atan2((double)0,-1) ; kappa = 0;

RotMat r =

1 -0 0

-0 -1 -1.22465e-016

0 1.22465e-016 -1

Quat q =

1 -0 0 -6.12323e-017

RotMat _r =

1 -0 0

0 -1 -1.22465e-016

0 1.22465e-016 -1

四元数k={{-1,0,0},{0,-1,0},{0,0,1}}

phi = 0 ; omiga = 0 ; kappa = atan2((double)0,-1);

RotMat r =

-1 -1.22465e-016 -0

1.22465e-016 -1 -0

0 -0 1

Quat q =

0 -0 1 -6.12323e-017

RotMat _r =

-1 -1.22465e-016 0

1.22465e-016 -1 -0

0 0 1

轴角到四元数的转换

假设旋转轴是(ax, ay, az)，记得必须是一个单位向量。

旋转角度是theta. （单位：弧度）

那么转换如下：

w = cos(theta / 2 )

x = ax * sin(theta / 2)

y = ay * sin(theta / 2)

z = az * sin(theta / 2 )

例如：

绕x轴旋转pi=>(1,0,0),pi，w=0,x=1,y=0,z=0。[w,(x,y,z)]=[0,(1,0,0)]=i

绕y轴旋转pi=>(0,1,0),pi，w=0,x=0,y=1,z=0。[w,(x,y,z)]=[0,(0,1,0)]=j

绕z轴旋转pi=>(0,0,1),pi，w=0,x=0,y=0,z=1。[w,(x,y,z)]=[0,(0,0,1)]=k

欧拉角到四元数的转换

如果你的欧拉角为(a, b, c)那么就可以形成三个独立的四元数，如下：

Qx=[ cos(a/2), (sin(a/2), 0, 0)]

Qy=[ cos(b/2), (0, sin(b/2), 0)]

Qz=[ cos(c/2), (0, 0, sin(c/2))]

最终的四元数是Qx * Qy * Qz的乘积的结果。
例如：

(0,0,0)=>Q_x=[1,(0,0,0)]，Q_y=[1,(0,0,0)]，Q_z=[1,(0,0,0)]，Q_xQ_yQ_z=1

(pi,0,0)=>Q_x=[0,(1,0,0)]，Q_y=[1,(0,0,0)]，Q_z=[1,(0,0,0)]，Q_xQ_yQ_z=Q_x=i

(0,pi,0)=>Q_x=[1,(0,0,0)]，Q_y=[0,(0,1,0)]，Q_z=[1,(0,0,0)]，Q_xQ_yQ_z=Q_y=j

(0,0,pi)=>Q_x=[1,(0,0,0)]，Q_y=[1,(0,0,0)]，Q_z=[0,(0,0,1)]，Q_xQ_yQ_z=Q_z=k

功能:四元变量的开方

格式:QuaternionSqrt(A)

a=[1 1 1 1]

a =

[ 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 ]

QuaternionSqrt(a)

ans =

[ 1.22474487139159 0.40824829046386 0.40824829046386 0.40824829046386 ]

功能:四元变量的次方

格式:QuaternionPow(A,B)

a=[1.22474487139159 0.40824829046386 0.40824829046386 0.40824829046386]

a =

[ 1.22474487139159 0.40824829046386 0.40824829046386 0.40824829046386 ]

QuaternionPow(a,2)

ans =

[ 1.00000000000001 0.99999999999999 0.99999999999999 0.99999999999999 ]