Lu数值计算扩展动态库LuMath

Lu数值计算扩展动态库LuMath 1.5

函数及类别	说明
1 什么是LuMath
2 LuMath数组（矩阵）对象基本用法
3 LuMath常量
math::real_s	用在函数new[math::real_s,... ...]中，申请数组或矩阵。该对象的基本类型为reals，扩展类型为real_s。
math::pi	math::pi=3.141592653589793
math::zeros	用在函数math::subg中，表示返回全0矩阵。
math::ones	用在函数math::subg中，表示返回全1矩阵。
math::eye	用在函数math::subg中，表示返回单位矩阵。
math::rand	用在函数math::subg中，表示返回0～1之间均匀分布的随机数矩阵。
math::randn	用在函数math::subg中，表示返回均值为0，方差为1的标准正态分布的随机数矩阵。
4 LuMath运算符
加法 +	数组（矩阵）加法。运行错误：1：不可识别对象；2：数组不存在；3：数组（矩阵）不匹配；4：内存错误。
减法 -	数组（矩阵）减法。运行错误：1：不可识别对象；2：数组不存在；3：数组（矩阵）不匹配；4：内存错误。
乘法 *	矩阵乘法。运行错误：1：不可识别对象；2：矩阵不存在；3：矩阵不匹配；4：内存错误。
转置 '	矩阵转置。运行错误：1：不可识别对象；2：矩阵不存在；3：内存错误。
点乘 .*	数组点乘。运行错误：1：不可识别对象；2：数组不存在；3：数组不匹配；4：内存错误。
点左除 ./	数组点左除。运行错误：1：不可识别对象；2：数组不存在；3：数组不匹配；4：内存错误。
点乘方 .^	数组点乘方。运行错误：1：不可识别对象；2：数组不存在；3：数组不匹配；4：内存错误。
5 LuMath函数
常用函数
math::toreal[cv,... ...]	将复数或三维向量转换为实数。转换成功返回true，否则返回false。 math::toreal[c,&x1,&x2] //c是一个复数，x1和x2接收复数的实部和虚部 math::toreal[v,&x1,&x2,&x3] //v是一个三维向量，x1、x2和x3接收三维向量的三个值
重载函数
sin、cos、... ...	该库对以下函数进行了重载，支持real_s类型对象的相关运算：new、copy、len、oset、oget、o、sqrt、exp、ln、lg、sin、cos、tan、asin、acos、atan、sinh、cosh、tanh、abs、floor、ceil、atan2、fmod。
in	与Lu核心库中用于简单序列化对象的in函数用法完全相同，可迭代获取一维数组、矩阵（当作一维数组进行操作）等序列化对象的值。
数（组）函数
math::ra1[2,3.5,-5]	申请一维数组并初始化。运行错误：1：不可识别参数，只接受实数或整数参数。
math::SetArray	对数组或矩阵进行初始化。
math::outa	输出数组或矩阵。
math::Sum	数组元素求和。
sum	该函数是重载函数，执行与math::Sum完全相同的运算。
math::ndgrid	生成n维函数的辅助格点阵列。
copy	该函数是重载函数。copy(a)将生成对象a的副本；copy(a,b)将数组b复制到数组a，两数组必须等长，并将多维数组视为一维维数；copy(a,abegin,b,bbegin,k)将数组b从地址bbegin开始的k个元素复制到数组a的自地址abegin开始的单元。运行错误：1：参数太多或者太少；2：数组不存在；3：数组不匹配；4：内存错误。
oset	该函数是重载函数。设a是一个二维数组对象，则执行形如a[i,j]=2.0的赋值语句时将调用该函数。运行错误：1：参数太少；2：不可识别对象；3：数组不匹配；4：参数不符合要求。
oget	该函数是重载函数。设a是一个二维数组对象，则执行形如a[i,j]的语句时将获得该对象元素的值。运行错误：1：参数太少；2：不可识别对象；3：数组不匹配；4：元素地址错误。
math::arrayfun	以实数数组元素为自变量进行函数调用。实际上，很容易实现以任意元素（lu元素）为自变量进行函数调用的函数lufun。
math::reshape	改变数组（矩阵）形状。
矩阵函数
math::matrix	申请m×n矩阵并初始化。
math::subg	返回一个矩阵的子矩阵。
math::subs	设置一个矩阵的子矩阵。
math::row(a,b,n) math::row(a,b,x)	产生一个单行矩阵（行向量）。a、b、x为实数，n为整数。a和b是行向量的第一个和最后一个元素，元素总数是n。或者，a是第一个元素，后面的元素按x递增，最后一个元素不超过b。
math::col(a,b,n) math::col(a,b,x)	产生一个单列矩阵（列向量）。a、b、x为实数，n为整数。a和b是列向量的第一个和最后一个元素，元素总数是n。或者，a是第一个元素，后面的元素按x递增，最后一个元素不超过b。
math::linspace(a,b,n) math::linspace(a,b,x)	产生一个一维数组。a、b、x为实数，n为整数。a和b是数组的第一个和最后一个元素，元素总数是n。或者，a是第一个元素，后面的元素按x递增，最后一个元素不超过b。
math::zeros(a) math::zeros(i,j)	将矩阵a重新初始化为全0矩阵，或者生成一个i×j全0矩阵。运行错误：1：数组不存在；2：不是二维数组（矩阵）；3：内存错误；4：参数太多或者太少。
math::ones(a) math::ones(i,j)	将矩阵a重新初始化为全1矩阵，或者生成一个i×j全1矩阵。运行错误：1：数组不存在；2：不是二维数组（矩阵）；3：内存错误；4：参数太多或者太少。
math::eye(a) math::eye(i,j)	将矩阵a重新初始化为单位矩阵，或者生成一个i×j单位矩阵。运行错误：1：数组不存在；2：不是二维数组（矩阵）；3：内存错误；4：参数太多或者太少。
math::rand(a) math::rand(i,j)	将矩阵a重新初始化为0～1之间均匀分布的随机数矩阵，或者生成一个i×j矩阵，并用0～1之间均匀分布的随机数初始化。运行错误：1：数组不存在；2：不是二维数组（矩阵）；3：内存错误；4：参数太多或者太少。
math::randn(a) math::randn(i,j)	将矩阵a重新初始化为均值为0，方差为1的标准正态分布的随机数矩阵，或者生成一个i×j矩阵，并用均值为0，方差为1的标准正态分布的随机数初始化。运行错误：1：数组不存在；2：不是二维数组（矩阵）；3：内存错误；4：参数太多或者太少。
math::inv(a)	矩阵求逆。运行错误：1：不可识别对象；2：不是n×n矩阵；3：内存错误；4：矩阵奇异。
解方程函数
math::netn	求非线性方程组一组实根的拟牛顿法。
math::btFindRoot	对分法求非线性方程的实根。
math::pqrt	求非线性方程一个实根的连分式解法。
积分函数
math::gsl_qng	QNG 非自适应gauss-krondrod积分。
math::gsl_qag	QAG 自适应积分。
math::gsl_qags	QAGS 带奇异值的自适应积分。
math::gsl_qagp	QAGP 考虑已知奇异点的自适应积分。
math::gsl_qagi	QAGI积分区间为（-∞，∞）的自适应积分。
math::gsl_qagiu	QAGIU积分区间为[a,∞）的自适应积分。
math::gsl_qagil	QAGI积分区间为（-∞，b]的自适应积分。
math::gsl_qawc	QAWC 柯西主值(Cauchy Principal Value)自适应积分
微分方程函数
math::gsl_ode	进化算法求解微分方程。

1 什么是LuMath [返回页首]

LuMath64.dll是一个Lu数值计算扩展动态库，该库以线性代数特别是矩阵运算为基础。

LuMath 1.5 由数学库GSL提供支持（仅此版本由GSL提供支持），下载源代码：http://www.forcal.net/xiazai/lu2/lugslmath64.rar

在LuMath中的函数是通过二级函数命名空间“math”输出的，所有函数均具有类似“math::matrix(...)”的格式。使用!!!using("math");可简化LuMath中的函数访问。

LuMath库的数组是C格式的，元素序号是基于0的。

LuMath库函数具有内存消耗低、执行效率高、代码简洁、实用性强的特点。

LuMath库可用于开发极致性能的应用程序，是熟悉C/C++、Fortran的数学爱好者的极佳的练手工具，同时也期望对一般的数值计算用户提供越来越多的方便。

限于作者水平，期待与朋友们共同完善LuMath！如果您有什么好的算法，任何改进的意见或建议，请与作者联系。

本库仅是抛砖引玉，任何喜欢LuMath的个人、团队或商业公司可基于此库开发商业程序。

2 LuMath数组（矩阵）对象基本用法 [返回页首]

借助Lu编译器的支持，LuMath中的数组（矩阵）对象有以下简洁的用法。

设A、B是LuMath数组（矩阵）对象：

A.=B //将对象B赋值给对象A

A.=A+B //将对象A+B赋值给对象A

A[2]=A[2]+1 //将一维数组A第2个单元的值加1

A[2,3,5]=A[3,1,2]+1 //将三维数组A的第[2,3,5]个单元的值赋值为A[3,1,2]+1

A.*2.0 //数组（矩阵）A的每一个元素乘以2.0

若A是一个LuMath矩阵：

A(type) //返回与矩阵A同样大小的矩阵，并根据type的值进行初始化。type的可选值为zeros（全0矩阵）、ones（全1矩阵）、eye（单位矩阵）、rand（0～1之间均匀分布的随机数）、randn（均值为0，方差为1的标准正态分布的随机数）。

A(all:j) //取矩阵A第j列所有元素。all是一个常量，在这里表示取所有行

A(i:all) //取矩阵A第i行所有元素。all是一个常量，在这里表示取所有列

A(i,i+m:all) //取矩阵A第i～i+m行所有元素。all是一个常量，在这里表示取所有列

A(all:j,j+m) //取矩阵A第j～j+m列所有元素。all是一个常量，在这里表示取所有行

A(i,i+m:j,j+n) //取矩阵A第i～i+m行，并在第j～j+n列的所有元素

A(all:j)=B //设置矩阵A第j列所有元素，B是一个数组。all是一个常量，在这里表示取所有行

A(i:all)=B //设置矩阵A第i行所有元素，B是一个数组。all是一个常量，在这里表示取所有列

A(i,i+m:all)=B //设置矩阵A第i～i+m行所有元素，B是一个矩阵。all是一个常量，在这里表示取所有列

A(all:j,j+m)=B //设置矩阵A第j～j+m列所有元素，B是一个矩阵。all是一个常量，在这里表示取所有行

A(i,i+m:j,j+n)=B //设置矩阵A第i～i+m行，并在第j～j+n列的所有元素，B是一个矩阵。

[例子]

!!!using["math"];
main(:t0,k,i,a,b)=
{
    t0=clock(),
    k=zeros(5,5),         //生成5×5矩阵k，初始化为0
    i=0,(++i<=1000).while{//循环计算1000次
        a=rand(5,7), b=rand(7,5), //生成5×7矩阵a，7×5矩阵b，用0～1之间的随机数初始化
        k.=k+a*b+a(0,4:1,5)*b(1,5:0,4)+a(all:6)*b(3:all)
    },
    k.outa(),             //输出矩阵k
    [clock()-t0]/1000.    //得到计算时间，秒
};

5 LuMath函数 [返回页首]

[返回页首] math::SetArray(a : x1,x2,... ...)：对数组或矩阵进行初始化

用法1：math::SetArray(a : x1,x2,... ...)：对数组或矩阵a进行初始化，从地址0开始存放数据x1,x2,... ...，初始化数据可以少于或多于数组大小，多余的数据被忽略。

用法2：math::SetArray(a : "1 -2 3 2.3 ... ...")：对数组或矩阵a进行初始化，从地址0开始存放数据，数据存放在字符串"1 -2 3 2.3 ... ..."中，初始化数据可以少于或多于数组大小，多余的数据被忽略。

运行错误：1：非法的数组对象；2：非法的Lu字符串及其他参数；3：非法的数字字符串。

[返回页首] math::outa(a, precision, rl, index)：输出数组或矩阵

    a：数组或矩阵，按列表方式输出。
    precision：有效数字位数（2～18），缺省是6位有效数字。
    rl：左对齐或右对齐输出，非0表示右对齐。缺省是右对齐。
    index：index非0时输出数组元素下标。缺省是不输出数组元素下标。

运行错误：1：参数太多或者太少；2：参数不符合要求；3：数组不存在；3：内存错误。

[返回页首] math::Sum(a,i)：数组元素求和

设有m维数组（含矩阵）：a(n1,n2,... ...,nm)，则Sum(a,i)对第i维求和，返回一个m-1维数组。若a是一维数组、1×k矩阵、k×1矩阵、i<1或者i>m，则Sum函数返回所有数组元素的和。Sum(a)相当于Sum(a,m)。

运行错误：1：参数太多或者太少；2：数组不存在；3：数组维数太多；4：内存错误。

[返回页首] math::ndgrid(x1,x2,...,xn : &y1,&y2,...,yn)：生成n维函数的辅助格点阵列

    x1,x2,...,xn：均是一维数组。
    &y1,&y2,...,yn：返回n维函数的辅助格点阵列。
    返回值：无意义。

运行错误：1：参数太多或者太少；2：数组不存在；3：不是一维数组；4：内存错误。

例子：

!!!using["math"];
mvar:
ndgrid[linspace(0.0,1.0,0.0011),linspace(1.0,2.0,0.0009),&x,&y],
Sum[cos(1.0-sin((1.2).*x.^(y./2.0)+cos(1.0-sin((1.2).*y.^(x./2.0))))),0];

相当于：

f(x,y)=cos(1.0-sin(1.2*x^(y/2.0)+cos(1-sin(1.2*y^(x/2.0)))));
sum[@f:0.0,1.0,0.0011:1.0,2.0,0.0009];

[返回页首] math::arrayfun(hFor:x1,x2,...,xn : xx)：以实数数组元素为自变量进行函数调用。

hFor为函数句柄；x1,x2,...,xn : xx为等长的实数数组（忽略维数）；xx用于存放结果，可缺省。若缺省xx，返回一个实数数组存放结果，否则该函数返回xx。该函数以数组x1,x2,...,xn的第i个元素为自变量，用函数hFor进行计算，结果存入xx的第i个单元。

运行错误：1：至少需要2个参数；2：指定的表达式不存在；3：数组个数与自变量个数不相等；4：不是实数数组；5：数组长度不匹配；6：内存错误；7：函数应返回一个实数。

例子：

!!!using["math"];
f(x,y)=x+y;
arrayfun[@f,linspace(1.,5.,5),linspace(1.,5.,5)].o[];

[返回页首] math::reshape(a,... ...)：改变数组（矩阵）形状

    格式1：math::reshape(a)：将矩阵重置为一维数组。
    格式2：math::reshape(a,i,j)：将矩阵a（或一维数组）重置为i×j矩阵，保持总元素不变。
    格式3：math::reshape(a,i,j,k,...)：改变多维数组的维数大小，保持维数和总元素不变。

运行错误代码：1：参数太多或者太少；2：实数数组不存在；3：数组维数小于1；4：数组不匹配。

[返回页首] math::matrix(m,n:x11,x12,... ...)：或者math::matrix(str)：申请m×n矩阵并初始化。

m,n：矩阵的行和列。初始化数据 x11,x12,...,x1n,x21,x22,...,x2n,... ...,xm1,xm2,...,xmn 可以缺省部分及全部，多余的数据被忽略。
str：字符串，包含矩阵数据，数据之间可以用空格或逗号分隔。此时，将自动确定矩阵的行和列。例如以下是一个2×3矩阵字符串：

"1 1.2 3,
2.2 3, 5"

运行错误：1：参数不符合要求；2：字符串矩阵每行中数据个数不同；3：字符串中数据非法。

[返回页首] math::subg(a, ... ...)：返回一个矩阵的子矩阵

    subg(a)：返回矩阵a的副本。
    subg(a:type)：返回与矩阵a同样大小的矩阵，并根据type的值进行初始化。type的可选值为zeros（全0矩阵）、ones（全1矩阵）、eye（单位矩阵）、rand（0～1之间均匀分布的随机数）、randn（均值为0，方差为1的标准正态分布的随机数）。
    subg(a:all:j)：取矩阵a第j列所有元素。all是一个常量，在这里表示取所有行。
    subg(a:i:all)：取矩阵a第i行所有元素。all是一个常量，在这里表示取所有列。
    subg(a:i,i+m:all)：取矩阵a第i～i+m行所有元素。all是一个常量，在这里表示取所有列。
    subg(a:all:j,j+m)：取矩阵a第j～j+m列所有元素。all是一个常量，在这里表示取所有行。
    subg(a:i,i+m:j,j+n)：取矩阵a第i～i+m行，并在第j～j+n列的所有元素。

运行错误：1：参数太多或者太少；2：数组不存在；3：不是二维数组（矩阵）；4：不可识别的类型标识type；5：无法定位矩阵的行或列；6：内存错误。

[例子]

!!!using["math"];
main(:a)=
    a=matrix[6,5].rand(),
    a.outa(),
    a.subg(all:3).outa(),
    a.subg(3:all).outa(),
    a.subg(3,5:2,3).outa();

[返回页首] math::subs(a, ... ...,b)：设置一个矩阵的子矩阵

    subs(a:b)：将矩阵b赋值给矩阵a，a与b应具有相同的结构。
    subs(a:all:j,b)：对矩阵a第j列所有元素赋值，b是一个等长的一维数组。all是一个常量，在这里表示取所有行。
    subs(a:i:all,b)：对矩阵a第i行所有元素赋值，b是一个等长的一维数组。all是一个常量，在这里表示取所有列。
    subs(a:i,i+m:all,b)：对矩阵a第i～i+m行所有元素赋值，b是一个相同大小的矩阵。all是一个常量，在这里表示取所有列。
    subs(a:all:j,j+m,b)：对矩阵a第j～j+m列所有元素赋值，b是一个相同大小的矩阵。all是一个常量，在这里表示取所有行。
    subs(a:i,i+m:j,j+n,b)：对矩阵a第i～i+m行，并在第j～j+n列的所有元素赋值，b是一个相同大小的矩阵。
    返回值：返回a。

运行错误：1：参数太多或者太少；2：数组不存在；3：不是二维数组（矩阵）；4：数组不匹配；5：无法定位矩阵的行或列。

[返回页首] math::netn(f,x,eps,t,h,k)：求非线性方程组一组实根的拟牛顿法

f：自定义二元或2n(n>1)元函数句柄，但必须由二级函数HFor获得该句柄。由用户自编。

如果f是二元函数，则两个参数为等长的一维数组：

f(x,y)=    //二元函数，x[i]为自变量，y[i]为方程右端点函数值
{
    y[0]=f1(x[0],x[1],...,x[n-1]),
    y[1]=f2(x[0],x[1],...,x[n-1]),
    ... ...
    y[n-1]=fn(x[0],x[1],...,x[n-1])
};

或者：

f(x1,x2,...,xn,y1,y2,...,yn)=    //2n元函数，xi为自变量，yi为方程右端点函数值
{
    y1=f1(x1,x2,...,xn),
    y2=f2(x1,x2,...,xn),
    ... ...
    yn=fn(x1,x2,...,xn)
};

x：双精度实型一维数组，长度不小于n，存放一组初值，返回方程组的一组实根。
eps：控制精度（eps>0）。可缺省该参数，缺省值eps=1e-6。
t：控制h大小的变量，0<t<1。可缺省该参数，缺省值h=0.1。
h：增量初值，在本函数中将被破坏。可缺省该参数，缺省值h=0.1。
k：允许的最大迭代次数。可缺省该参数，缺省值k=100。

返回值：返回值为实际迭代次数。若返回值为-1表示代数方程奇异，若返回值为-2表示β=0，这两种情况可放宽精度要求、改变各个初值或改变各个方程顺序再试；若返回值等于0说明迭代了k次仍未满足精度要求，程序工作失败。

运行错误代码：1：参数太多或者太少；2：指定的表达式不存在；3：方程f参数个数不符合要求；4：实数数组不存在；5：数组元素太少；6：精度、控制变量、迭代次数等不符合要求；7：内存错误；8：用户删除了自定义二元函数中的工作数组。

[例子]

f(x1,x2,x3,y1,y2,y3)=      //函数定义
{
    y1=x1*x1+x2*x2+x3*x3-1.0,
    y2=2.0*x1*x1+x2*x2-4.0*x3,
    y3=3.0*x1*x1-4.0*x2+x3*x3
};
!!!using["math"];
main(:x,i)={
    x=new{real_s,data : 1.0,1.0,1.0},   //设置初值为1.0,1.0,1.0
    i=netn[@f,x],          //拟牛顿法解方程
    x.outa(),             //输出结果
    i                      //返回迭代次数
};

或者：

f(x,y)=                   //函数定义
{
    y[0]=x[0]*x[0]+x[1]*x[1]+x[2]*x[2]-1.0,
    y[1]=2.0*x[0]*x[0]+x[1]*x[1]-4.0*x[2],
    y[2]=3.0*x[0]*x[0]-4.0*x[1]+x[2]*x[2]
};
!!!using["math"];
main(:x,i)={
    x=new{real_s,data : 1.0,1.0,1.0},   //设置初值为1.0,1.0,1.0
    i=netn[@f,x],          //拟牛顿法解方程
    x.outa(),             //输出结果
    i                      //返回迭代次数
};

[返回页首] math::btFindRoot(f,a,b,h,k,eps)：对分法求非线性方程的实根

f：自定义一元函数句柄。函数自变量不能重新赋值。
a,b：求根区间的左端点和右端点（b>a）。
h：搜索求根时采用的步长（h>0）。
k: 可能的实根个数。可缺省，缺省值为10。
eps：控制精度（eps>0）。可缺省，缺省值为1e-6。

返回值：返回一个数组（存放搜索到的实根），或者返回nil。若返回数组的长度为k，则有可能在求根区间[a,b]内的实根未搜索完。

运行错误代码：1：参数太多或者太少；2：指定的函数不存在；3：不是一元方程；4：自变量不能重新赋值；5：参数设置不符合要求；6：内存错误；7：指定的函数应返回一个实数。

例子：

f(x)=x^6-5*x^5+3*x^4+x^3-7*x^2+7*x-20; //函数定义
math::btFindRoot[@f,-2.,5.,0.2].o[];

[返回页首] math::qrrt(f,&x,eps)：求非线性方程一个实根的连分式解法

f：自定义一元函数句柄。
x：存放迭代初值，返回时存放迭代终值。
eps：控制精度（eps>0）。

返回值：返回值为实际迭代次数。若等于10，表示已迭代10次，但未满足精度要求，返回值仅作参考；若返回值小于10表示正常返回。

运行错误代码：1：指定的表达式不存在；2：不是一元方程；3：函数自变量重新赋值；4：参数不符合要求；5：函数应返回一个实数。

例子：

f(x)=x^3-x^2-1; //函数定义
main(:x,i)=
{
    x=1.0,
    i=math::pqrt[@f,&x,1e-6],
    o{i, " ", x, "\r\n"}
};

[返回页首] math::gsl_qng(f,a,b,epsabs,epsrel,&abserr,&errcode)：QNG非自适应gauss-krondrod积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。
a,b：积分上下限，均为实数。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

运行错误：1：嵌套太多调用；2：参数太多或者太少；3：指定的表达式不存在；4：不是一元函数；5：参数不符合要求；6：自定义一元函数应返回一个实数。

例子：

f(x)=ln(x)/sqrt(x);
math::gsl_qng[@f,0.1,1.]; //得到积分值
(:abserr)=math::gsl_qng[@f,0.1,1.,1e-6,1e-6,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qag(f,a,b,epsabs,epsrel,limit,key,&abserr,&errcode)：QAG自适应积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。
a,b：积分上下限，均为实数。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
limit：子区间划分上限，取值范围为5～100。可缺省，默认为10。
key：高斯积分点选项。可缺省，默认为4（gsl_gauss41）。
    gsl_gauss15 = 1 /* 15 point Gauss-Kronrod rule */
    gsl_gauss21 = 2 /* 21 point Gauss-Kronrod rule */
    gsl_gauss31 = 3 /* 31 point Gauss-Kronrod rule */
    gsl_gauss41 = 4 /* 41 point Gauss-Kronrod rule */
    gsl_gauss51 = 5 /* 51 point Gauss-Kronrod rule */
    gsl_gauss61 = 6 /* 61 point Gauss-Kronrod rule */
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

f(x)=ln(x)/sqrt(x);
math::gsl_qag[@f,0.1,1.]; //得到积分值
(:abserr)=math::gsl_qag[@f,0.1,1.,1e-6,1e-6,10,4,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qags(f,a,b,epsabs,epsrel,limit,&abserr,&errcode)：QAGS 带奇异值的自适应积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。
a,b：积分上下限，均为实数。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
limit：子区间划分上限，取值范围为5～100。可缺省，默认为10。
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

f(x)=ln(x)/sqrt(x);
math::gsl_qags[@f,0.1,1.]; //得到积分值
(:abserr)=math::gsl_qags[@f,0.1,1.,1e-6,1e-6,10,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qagp(f,math::ra1[a, x1, x2, ... , b],epsabs,epsrel,limit,&abserr,&errcode)：QAGP 考虑已知奇异点的自适应积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。
ra1[a, x1, x2, ... , b]：积分上下限a、b及已知奇异点x1, x2, ...。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
limit：子区间划分上限，取值范围为5～100。可缺省，默认为10。
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

f(x)=x^3*ln{abs[(x^2-1)*(x^2-2)]}; //函数在1和sqrt(2.0)处奇异
math::gsl_qagp[@f,math::ra1[0,1,sqrt(2.0),3],1e-6,1e-6,30]; //得到积分值
(:abserr)=math::gsl_qagp[@f,math::ra1[0,1,sqrt(2.0),3],1e-6,1e-6,30,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qagi(f,epsabs,epsrel,limit,&abserr,&errcode)：QAGI积分区间为（-∞，∞）的自适应积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
limit：子区间划分上限，取值范围为5～100。可缺省，默认为10。
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

f(x)=1/(x*x+1);
math::gsl_qagi[@f]; //得到积分值
(:abserr)=math::gsl_qagi[@f,1e-6,1e-6,10,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qagiu(f,a,epsabs,epsrel,limit,&abserr,&errcode)：QAGIU积分区间为[a,∞)的自适应积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。
a：积分下限，实数。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
limit：子区间划分上限，取值范围为5～100。可缺省，默认为10。
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

f(x)=1/(x*x+1);
math::gsl_qagiu[@f,1.]; //得到积分值
(:abserr)=math::gsl_qagiu[@f,1.,1e-6,1e-6,10,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qagil(f,a,epsabs,epsrel,limit,&abserr,&errcode)：QAGI积分区间为(-∞，b]的自适应积分，嵌套调用层数最多为6层。

返回值：积分值。

f(x)=l/(x*x+1);
math::gsl_qagil[@f,1.]; //得到积分值
(:abserr)=math::gsl_qagil[@f,1.,1e-6,1e-6,10,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_qawc(f,a,b,c,epsabs,epsrel,limit,&abserr,&errcode)：QAWC 柯西主值(Cauchy Principal Value)自适应积分，嵌套调用层数最多为6层。

f：自定义一元函数句柄。函数自变量不能重新赋值。该函数的形式为 F(x)=f(x)/(x-c)，仅取f(x)部分。
a,b,c：积分上下限a、b及带有一个奇异点c的柯西主值，均为实数。
epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
limit：子区间划分上限，取值范围为5～100。可缺省，默认为10。
abserr: 返回绝对误差值估计值。可缺省。
errcode：返回运行错误代码，返回0表示工作正常。可缺省。

返回值：积分值。

f(x)=exp(x);
math::gsl_qawc[@f,1.0,2.0,3.0]; //得到积分值
(:abserr)=math::gsl_qawc[@f,1.0,2.0,3.0,1e-6,1e-6,10,&abserr],abserr; //得到积分值的绝对误差

[返回页首] math::gsl_ode(f, jac, params, math::ra1[t0, t1, ..., tk], math::ra1[Y1,Y2,...,Yn], epsabs, epsrel, type, step, max, errtype)：进化算法求解微分方程，嵌套调用层数最多为6层。

f：用户定义的函数，用于计算微分方程组中各方程右端函数值，由用户自编。该表达式有n*n+2*n+2个参数，第一个参数为自变量，随后n个参数为函数值，随后n个参数为右端函数值（即微分方程的值），最后一个参数为gsl_ode的输入参数params。n为微分方程组中方程的个数，也是未知函数的个数。该函数计算成功时，应该返回0。该形式如下：

f(t,y1,y2,y3,d1,d2,d3,params)=
{
    d1 = y2,    //d(i) = f(t,yi,params)
    d2 = -y1,
    d3 = -y3,

    0 //计算成功返回0
};

jac：用户定义的函数，用于计算Jacobian矩阵，由用户自编。该表达式有2*n+2个参数，第一个参数为自变量，随后n个参数为函数值，随后n*n个参数为Jacobian 矩阵, 存储方式为J(i,j) = dfdy[i * n+ j]，随后n个参数为∂fi/∂t，最后一个参数为gsl_ode的输入参数params。n为微分方程组中方程的个数，也是未知函数的个数。对于一些简单的算法并不需要 Jacobian 矩阵，此时可以将该参数设为 nil，但是大多数高效的算法都要用到 Jacobian 矩阵，因此在能提供 Jacobian 矩阵时应该尽量提供。该函数计算成功时，应该返回0。该形式如下：

jac(t,y1,y2,y3,J11,J12,J12,J21,J22,J23,J31,J32,J33,d1,d2,d3,params)=
{
    J11 = ... ...,    //J(i,j) = f(t,yi,params)
    J12 = ... ...,
    J13 = ... ...,
    J21 = ... ...,
    J22 = ... ...,
    J23 = ... ...,
    J31 = ... ...,
    J32 = ... ...,
    J33 = ... ...,

    d1 = ... ...,    //d(i) = f(t,yi,params)
    d2 = ... ...,
    d3 = ... ...,

    0 //计算成功返回0
};

    params：输入参数，供自定义函数 f 和 jac 使用。如果 params = nil，将传入下一个参数 math::ra1[t0, t1, ..., tk] 当前数组元素的序号（整数），此序号是基于0的；例如序号为2时，表示该函数正在t2和t3之间求解。
    math::ra1[t0, t1, ..., tk]：数组，存放独立变量t。输入初始时间为t0，用户希望得到时间t1, ..., tk时的解。
    math::ra1[Y1,Y2,...,Yn]：存放n个未知函数在起始点处的函数值Y(n)=Yn(t)。
    epsabs,epsrel：绝对误差上限，相对误差上限，均为不小于0的实数。可缺省，缺省值均为1e-6。
    type：求解函数类型。可缺省该参数，缺省值1（gsl_rk4）。
        gsl_rk2   0 //二阶 Runge-Kutta 法
        gsl_rk4   1 //四阶经典 Runge-Kutta 法
        gsl_rkf45 2 //Runge-Kutta-Fehlberg 法
        gsl_rkck   3 //Runge-Kutta Cash-Karp
        gsl_rk8pd 4 //Runge-Kutta Prince-Dormand
        gsl_rk2imp 5 //隐式的 Runge-Kutta
        gsl_rk2simp 6 //未用
        gsl_rk4imp 7 //Implicit 4th order Runge-Kutta at Gaussian points
        gsl_bsimp 8 //隐式的 Bulirsch-Stoer method of Bader and Deuflha （需要 Jacobian 矩阵）--未用
        gsl_gear1 9 //M=1 implicit Gear method
        gsl_gear2 10 //M=2 implicit Gear method
    step：初始步长，该步长会被自动调整。可缺省该参数，缺省值为1e-6。
    max：允许最大迭代次数。可缺省该参数，缺省值为10000。
    errtype：达到最大迭代次数max时返回的错误类型，该参数为逻辑值。errtype=true时返回运行错误，错误类型为1；否则返回运行警告，警告类型为-2。可缺省该参数，缺省值为false。

返回值：(1+k)×(1+n)二维数组。存放(1+k)组结果（包含初始值），每一组结果即：t,Y1,Y2,...,Yn。

运行错误：1：不允许嵌套调用；2：参数太多或者太少；3：指定的表达式不存在；4：函数参数个数不符合要求；5：参数不符合要求；6：内存错误；7：函数应返回一个实数。

    [例子1] 设一阶微分方程组及初值为：
        r'=2r-2rf, r(0)=1
        f'=-f+rf, f(0)=3
    计算t=1,2,...,10时的r、f的值。

    程序如下：

!!!using["math"];
f(t,r,f,dr,df, p)={dr=2*r-2*r*f, df=-f+r*f, 0}; //函数定义
gsl_ode[@f,nil,0.0,ra1[0,1,2,3,4,5,6,7,8,9,10],ra1[1,3]].o[]; //不使用 Jacobian 矩阵

    [例子2] 考虑二阶的非线性 VanderPol 方程：

    x′′(t)+μx′(t)(x(t)^2−1)+x(t)=0

    这个方程可以被化简为如下的常微分方程组：

    y1' = -y0 + μ *y1* (1 - y0^2)
    y0' = y1

    Jacobian 矩阵为：

    0                         1
    −1−2*μ*y0*y1   μ*(1−y0^2)

设μ=10，已知t=0时的初值为[1, 0]，计算t=10,20,30,40,50,60,70,80,90,100时的y0、y1的值。

程序如下：

!!!using["math"];
f(t,y0,y1,f0,f1, p)=
{
    f0 = y1,
    f1 = -y0 - p*y1*(y0*y0 - 1),
    0
};
jac(t, y0,y1, j00,j01,j10,j11, f0,f1, p)=
{
    j00 = 0.0,
    j01 = 1.0,
    j10 = -2.0*p*y0*y1 - 1.0,
    j11 = -p*(y0*y0 - 1.0),

    f0 = 0.0,
    f1 = 0.0,

    0
};
gsl_ode[@f,@jac,10.0,ra1[0,10,20,30,40,50,60,70,80,90,100],ra1[1,0]].o[];