欢迎访问 Lu程序设计

Lu参数优化动态库LuOpt V2.0

目  录

1 什么是LuOpt [返回页首]

    LuOpt64.dll是一个Lu参数优化动态库，主要包含一些参数优化函数和解方程函数。

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

    免费注册说明：为了演示设计商业性Lu扩展动态库的可行性，LuOpt设计成须注册后使用其全部功能。如果没有注册，当优化参数多于5个时，Opt函数仅返回最小值，不能返回优化参数。目前免费注册，用函数luopt::OptPassword获取机器码后，通过E-mail发送到forcal@sina.com获取注册码 （已取消注册，可直接使用LuOpt的全部功能）。

    函数OptPassword的格式（返回值及参数pw1,pw1,... ...,pwn将返回多个机器码）：luopt::OptPassword(&pw1,&pw1,... ...,&pwn)

    函数OptPassword的用法1：luopt::OptPassword();

    函数OptPassword的用法2：(:pw)=luopt::OptPassword(&pw),pw;

    函数OptPassword的用法3：(:pw)=luopt::OptPassword(0,&pw),pw;

    获取注册码后，在加载LuOpt64.dll时提供注册码即可。一般格式为（注册码在冒号后）："dll\LuOpt64.dll:604508320"

2 Lu参数优化函数 [返回页首]

[返回页首] luopt::Opt(a, f, luopt::optconfun,r, luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax, luopt::optautorange,autorange, luopt::optmin0,min0, luopt::optmax,max, luopt::optbuf,&buf, luopt::optmode,mode, luopt::optdeep,deep, luopt::optbreadth,breadth, luopt::optkeepmins,keepmins, luopt::optfast,fast, luopt::optwaynull, luopt::optwaysim, luopt::optwaysimdeep, luopt::optwayconfra, luopt::optwaygrad,grad, luopt::optwaygrow, luopt::optwaycom, luopt::opteps,eps, luopt::optsameeps,sameeps, luopt::optsameabs,sameabs, luopt::optout,out, luopt::optshow,show, luopt::optnum,num, luopt::optarray,a, luopt::optarraybegin,&k, luopt::optarraystarter,arraystarter, luopt::optexpand,expand, luopt::optcontract,contract, luopt::optsimratio,simratio,  luopt::optstarter,&x1, &x2, ... ... , &xn, &min)：求无约束或约束条件下的n维极小值

a： 一维实数数组，数组长度n即拟合参数的个数；该数组将作为目标函数f的输入参数，目标函数f为一元函数；该参数可以缺省，缺省时目标函数可为n元函数，n元函数的自变量即n个拟合参数。

f：目标函数句柄，不可缺省，该函数由用户自编。如果f是Opt的第一个参数，f为自定义n元函数句柄，用于计算目标函数f(x1,x2,...,xn)的值，格式如下：

f(x1,x2,...,xn)=
{
    g(x1,x2,...,xn)
};

    如果f是第二个参数，f为自定义 一元函数句柄，用于计算目标函数f(a)的值，其中a为一维实数数组，即Opt的第一个参数，格式如下：

f(x)=
{
    g[x(1),x(2),...,x(n)]
};

luopt::optconfun,r： 约束函数句柄，可以缺省，该函数由用户自编。约束函数r的参数个数及格式与目标函数相同，如果定义了参数a，约束函数为一元函数；如果没有定义参数a，约束函数为n元函数。约束函数返回非0值表示满足约束条件，若不满足约束条件，应返回0、false或nil。

luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax：指定搜索区间 。搜索区间越小，越容易获得最优值。若缺省该参数，则所有变量区间均为[-1e20～1e20]。
luopt::optautorange,autorange：autorange!=0时将在指定搜索区间内自动调节合适的搜索区间，否则不进行调节。默认是自动调节搜索区间的。
luopt::optmin0,min0：解的可能的极小值。可以缺省该参数，缺省值为±1e-20。
luopt::optmax,max：max>=10，控制搜索精度，max越大精度越高，但时间越长。可以缺省该参数，缺省值为100。
luopt::optbuf,buf：buf>=1，控制缓冲区大小，缓冲区并非越大越好。可以缺省该参数，缺省值为5。
luopt::optmode,mode：mode>=0，控制工作模式。mode越大越容易搜索到最优值。可以缺省该参数，缺省值为16。
luopt::optdeep,deep：deep>=0，控制搜索深度。deep越大越容易搜索到最优值。可以缺省该参数，缺省值为20。 当deep为偶数或奇数时使用不同的搜索策略；当deep为奇数时搜索相似点。
luopt::optbreadth,breadth：breadth>=0，控制搜索广度以替代默认的搜索方法。breadth越大搜索范围越大。 通常breadth取10～20。可以缺省该参数，breadth=0也相当于缺省该参数。
luopt::optkeepmins,keepmins：keepmins>=0，控制下一轮搜索时保留的极小值数目。可以缺省该参数，keepmins=0也相当于缺省该参数 ，此时保留所有极小值。
luopt::optfast,fast：fast!=0，进行快速搜索。可以缺省该参数，缺省值为fast=0。
luopt::optwaynull：清除所有局部搜索方法。
luopt::optwaysim：使用单形调优法局部搜索，这是默认的局部搜索方法。
luopt::optwaylm：使用麦夸特法（Levenberg-Marquardt）局部搜索。
luopt::optwaylme：使用扩展式麦夸特法（Levenberg-Marquardt）局部搜索，本方法包含了optwaylm方法。
luopt::optwaynewton：使用牛顿法局部搜索。
luopt::optwaysimdeep：使用单形调优法局部深度搜索。
luopt::optwayconfra：使用连分式局部搜索方法。
luopt::optwaygrad,grad：使用梯度搜索方法。grad>=0，控制梯度搜索的深度。grad越大耗时越长。通常取grad=0。
luopt::optwaygrow：使用扩展式搜索方法。参数多时该方法较慢。
luopt::optwaycom：使用复形调优法局部搜索方法。用optconfun指定约束函数时自动启用该方法。
luopt::opteps,eps：控制精度要求，eps>0。可以缺省该参数，缺省值为1e-6。
luopt::optsameeps,sameeps：相对比较两组值是否相等的极小值。sameeps>0。可以缺省该参数，缺省值为1e-2。
luopt::optsameabs,sameabs：绝对比较两组值是否相等的极小值。sameeps>=0。可以缺省该参数，缺省值为1e-3。
luopt::optout,out：是否输出结果。out非0时输出。可以缺省该参数，默认输出结果。
luopt::optshow,show：是否通过窗口输出中间结果。show非0时输出。可以缺省该参数，默认输出中间结果。
luopt::optnum,num：设置最多输出的极小值的个数。可以缺省该参数，默认输出1个极小值。
luopt::optarray,a：提供一个数组，用于接受结果。一组结果包括所有的自变量参数及极小值。默认从数组开始存放数据，也可用参数optarraybegin指定开始存放数据的位置。可以缺省该参数。
luopt::optarraybegin,&k：指定数组中开始存放数据的位置，k返回下一个存放位置。可以缺省该参数，默认从地址0开始存放。
luopt::optarraystarter,arraystarter：arraystarter!=0将指定数组中包含一组初值。可以缺省该参数，缺省值为arraystarter=0。若同时指定了optstarter，optstarter中的初值被忽略。
luopt::optexpand,expand：单形调优法扩张系数。一般取1.2～2.0。可以缺省该参数，缺省值为1.6。
luopt::optcontract,contract：单形调优法收缩系数。一般取0～1。可以缺省该参数，缺省值为0.5。
luopt::optsimratio,simratio：单形调优法调整系数。一般取0～2之间的数，但不可取0和1这两个值。可以缺省该参数，缺省值为0.9。
luopt::optstarter,&x1, &x2, ... ... , &xn, &min：前n个参数提供一组初值（可为任意值），返回最优参数值 ，最后一个参数返回极小值。可以缺省该参数。

返回值：极小值个数。

说明1：luopt::optmax、luopt::optbuf、luopt::opteps等标识参数类型，次序是任意的，可缺省。

说明2：该函数使用随机算法搜索全局最小值，故解是不稳定的，应多次搜索甚至变换参数搜索，以最小者为最优。

说明3：该函数将自动打开一个窗口输出中间结果，用户可通过分析这些中间结果以决定是否继续求解。

运行错误：

1：指定的表达式不存在；
2：常量表达式或者自变量不能重新赋值；
3：内存错误；
4：参数不匹配；
5：不可识别描述符；
6：参数不符合要求；
7：约束函数有错误：参数不匹配、自变量重新赋值或者函数不存在；
8：自定义函数返回了非法值；
9：该函数不允许嵌套调用。

例子1：计算下列目标函数的极小值点与极小值。
    J=100*(x1-x0*x0)²+(1-x0)²

程序如下：

f(x0,x1)=100*(x1-x0*x0)^2+(1-x0)^2;
luopt::Opt[@f];

结果：

例子2：求下列隐函数z的最小值：

    z=sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}

    其中x范围[-1,7]，y范围[-2,2]。

程序如下：

!!!using("luopt");
f(x,y,z)=z+1e10*{z-sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}}^2; //构造目标函数，注意1e10
Opt[@f, optwaylm, optrange,-1.,7., -2.,2., -1e10,1e10];
//Opt[@f, optwaycom, optrange,-1.,7., -2.,2., -1e10,1e10];
//Opt[@f, optwaynewton, optrange,-1.,7., -2.,2., -1e10,1e10];    //也可以使用此语句
//Opt[@f, optwaysimdeep, optwayconfra, optrange,-1.,7., -2.,2., -1e10,1e10];    //也可以使用此语句
//Opt[@f, optwaygrow, optrange,-1.,7., -2.,2., -1e10,1e10];    //也可以使用此语句

结果：

2.898270235856681 -0.8573130005839527 -2.335408331097716e-002 -2.335408328621853e-002

如果不指定范围，程序如下：

!!!using("luopt");
f(x,y,z)=z+1e10*{z-sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}}^2; //构造目标函数，注意1e10
Opt[@f, optwaygrow, optwaygrad,0, optmode,200];

结果（需多次运行）：

2.898270235856681 -0.8573130005839527 -2.335408331097716e-002 -2.335408328621853e-002

例子3：求针状函数的全局最小值，并返回10个极小值：

    f(r)=-sin(r)/r+1

    其中:r=sqrt[(x-50)^2+(y-50)^2]+exp[1.0]。

程序如下：

!!!using("luopt");
f(x,y:t)=   //目标函数定义
{
    t=sqrt[(x-50)^2+(y-50)^2]+exp[1.0],
    -sin[t]/t-1
};
Opt[@f, optwaygrad,0, optnum,10];

结果（需多次运行）：

50.               50.               -1.151117991593894
45.04713670154168 50.73409357987649 -1.128374553525899
47.95358156651113 54.56960337517742 -1.128374553259392
48.61084727988142 45.18968474411513 -1.128374553040821
51.53888994360668 55.27956633578699 -1.11373603845343
50.0914003285785  38.652456006392   -1.070913459450462
55.57756844733133 40.1173971653779  -1.070913459450462
47.78264700546625 38.87080820770154 -1.070913459436176
43.59630987066741 40.6315049957481  -1.070913459402211
46.75204544703252 39.12675998582707 -1.070913459306036
10.

[返回页首] luopt::Opt1(a, f, luopt::optconfun,r, luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax, luopt::optmin0,min0, luopt::optmax,max, luopt::optbuf,&buf, luopt::optmode,mode, luopt::optdeep,deep, luopt::optwaynull, luopt::optwaysim, luopt::optwaysimdeep, luopt::optwayconfra, luopt::optwaygrad,grad, luopt::optwaygrow, luopt::optwaycom, luopt::opteps,eps, luopt::optsameeps,sameeps, luopt::optsameabs,sameabs, luopt::optout,out, luopt::optshow,show, luopt::optnum,num, luopt::optarray,a, luopt::optarraybegin,&k, luopt::optarraystarter,arraystarter, luopt::optexpand,expand, luopt::optcontract,contract, luopt::optsimratio,simratio,  luopt::optstarter,&x1, &x2, ... ... , &xn, &min)：求无约束或约束条件下的n维极小值

a： 一维实数数组，数组长度n即拟合参数的个数；该数组将作为目标函数f的输入参数，目标函数f为一元函数；该参数可以缺省，缺省时目标函数可为n元函数，n元函数的自变量即n个拟合参数。

f：目标函数句柄，不可缺省，该函数由用户自编。如果f是Opt1的第一个参数，f为自定义n元函数句柄，用于计算目标函数f(x1,x2,...,xn)的值，格式如下：

f(x1,x2,...,xn)=
{
    g(x1,x2,...,xn)
};

    如果f是第二个参数，f为自定义 一元函数句柄，用于计算目标函数f(a)的值，其中a为一维实数数组，即Opt的第一个参数，格式如下：

f(x)=
{
    g[x(1),x(2),...,x(n)]
};

luopt::optconfun,r： 约束函数句柄，可以缺省，该函数由用户自编。约束函数r的参数个数及格式与目标函数相同，如果定义了参数a，约束函数为一元函数；如果没有定义参数a，约束函数为n元函数。约束函数返回非0值表示满足约束条件，若不满足约束条件，应返回0、false或nil。

luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax：指定搜索区间 。搜索区间越小，越容易获得最优值。若缺省该参数，则所有变量区间均为[-1e20～1e20]。
luopt::optmin0,min0：解的可能的极小值。可以缺省该参数，缺省值为±1e-20。
luopt::optmax,max：max>=10，控制搜索精度，max越大精度越高，但时间越长。可以缺省该参数，缺省值为100。
luopt::optbuf,buf：buf>=1，控制缓冲区大小，缓冲区并非越大越好。可以缺省该参数，缺省值为3。
luopt::optmode,mode：mode>=0，控制工作模式。mode越大越容易搜索到最优值。可以缺省该参数，缺省值为6。当mode为偶数或奇数时使用不同的搜索策略；当mode为奇数时更重视在极小值点处搜索最优点。
luopt::optdeep,deep：deep>=0，控制搜索深度。deep越大越容易搜索到最优值。可以缺省该参数，缺省值为6。当deep为偶数或奇数时使用不同的搜索策略；当deep为奇数时搜索相似点。
luopt::optwaynull：清除所有局部搜索方法。
luopt::optwaysim：使用单形调优法局部搜索，这是默认的局部搜索方法。
luopt::optwaylm：使用麦夸特法（Levenberg-Marquardt）局部搜索。
luopt::optwaylme：使用扩展式麦夸特法（Levenberg-Marquardt）局部搜索，本方法包含了optwaylm方法。
luopt::optwaynewton：使用牛顿法局部搜索。
luopt::optwaysimdeep：使用单形调优法局部深度搜索。
luopt::optwayconfra：使用连分式局部搜索方法。
luopt::optwaygrad,grad：使用梯度搜索方法。grad>=0，控制梯度搜索的深度。grad越大耗时越长。通常取grad=0。
luopt::optwaygrow：使用扩展式搜索方法。参数多时该方法较慢。
luopt::optwaycom：使用复形调优法局部搜索方法。用optconfun指定约束函数时自动启用该方法。
luopt::opteps,eps：控制精度要求，eps>0。可以缺省该参数，缺省值为1e-6。
luopt::optsameeps,sameeps：相对比较两组值是否相等的极小值。sameeps>0。可以缺省该参数，缺省值为1e-2。
luopt::optsameabs,sameabs：绝对比较两组值是否相等的极小值。sameeps>=0。可以缺省该参数，缺省值为1e-3。
luopt::optout,out：是否输出结果。out非0时输出。可以缺省该参数，默认输出结果。
luopt::optshow,show：是否通过窗口输出中间结果。show非0时输出。可以缺省该参数，默认输出中间结果。
luopt::optnum,num：设置最多输出的极小值的个数。可以缺省该参数，默认输出1个极小值。
luopt::optarray,a：提供一个数组，用于接受结果。一组结果包括所有的自变量参数及极小值。默认从数组开始存放数据，也可用参数optarraybegin指定开始存放数据的位置。可以缺省该参数。
luopt::optarraybegin,&k：指定数组中开始存放数据的位置，k返回下一个存放位置。可以缺省该参数，默认从地址0开始存放。
luopt::optarraystarter,arraystarter：arraystarter!=0将指定数组中包含一组初值。可以缺省该参数，缺省值为arraystarter=0。若同时指定了optstarter，optstarter中的初值被忽略。
luopt::optexpand,expand：单形调优法扩张系数。一般取1.2～2.0。可以缺省该参数，缺省值为1.6。
luopt::optcontract,contract：单形调优法收缩系数。一般取0～1。可以缺省该参数，缺省值为0.5。
luopt::optsimratio,simratio：单形调优法调整系数。一般取0～2之间的数，但不可取0和1这两个值。可以缺省该参数，缺省值为0.9。
luopt::optstarter,&x1, &x2, ... ... , &xn, &min：前n个参数提供一组初值（可为任意值），返回最优参数值 ，最后一个参数返回极小值。可以缺省该参数。

返回值：极小值个数。

说明1：luopt::optmax、luopt::optbuf、luopt::opteps等标识参数类型，次序是任意的，可缺省。

说明2：该函数使用随机算法搜索全局最小值，故解是不稳定的，应多次搜索甚至变换参数搜索，以最小者为最优。

说明3：该函数将自动打开一个窗口输出中间结果，用户可通过分析这些中间结果以决定是否继续求解。

运行错误：

1：指定的表达式不存在；
2：常量表达式或者自变量不能重新赋值；
3：内存错误；
4：参数不匹配；
5：不可识别描述符；
6：参数不符合要求；
7：约束函数有错误：参数不匹配、自变量重新赋值或者函数不存在；
8：自定义函数返回了非法值；
9：该函数不允许嵌套调用。

例子1：计算下列目标函数的极小值点与极小值。
    J=100*(x1-x0*x0)²+(1-x0)²

程序如下：

f(x0,x1)=100*(x1-x0*x0)^2+(1-x0)^2;
luopt::Opt1[@f];

结果：

例子2：求下列隐函数z的最小值：

    z=sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}

    其中x范围[-1,7]，y范围[-2,2]。

程序如下：

!!!using("luopt");
f(x,y,z)=z+1e10*{z-sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}}^2; //构造目标函数，注意1e10
Opt1[@f, optwaycom, optrange,-1.,7., -2.,2., -1e10,1e10];
//Opt1[@f, optwaylm, optrange,-1.,7., -2.,2., -1e10,1e10];
//Opt1[@f, optwaynewton, optrange,-1.,7., -2.,2., -1e10,1e10];    //也可以使用此语句
//Opt1[@f, optmax,10, optmode,1, optbuf,1, optdeep,1, optwaynewton, optrange,-1.,7., -2.,2., -1e10,1e10];    //也可以使用此语句

结果：

2.898270235856681 -0.8573130005839527 -2.335408331097716e-002 -2.335408328621853e-002

如果不指定范围，程序如下：

!!!using("luopt");
f(x,y,z)=z+1e10*{z-sin[(z*x-0.5)^2+2*x*y*y-z/10]*exp{-[(x-0.5-exp(-y+z))^2+y*y-z/5+3]}}^2; //构造目标函数，注意1e10
Opt1[@f, optwaygrow, optwaygrad,0, optmode,50];

结果（需多次运行）：

2.898270235856681 -0.8573130005839527 -2.335408331097716e-002 -2.335408328621853e-002

例子3：求针状函数的全局最小值，并返回10个极小值：

    f(r)=-sin(r)/r+1

    其中:r=sqrt[(x-50)^2+(y-50)^2]+exp[1.0]。

程序如下：

!!!using("luopt");
f(x,y:t)=   //目标函数定义
{
    t=sqrt[(x-50)^2+(y-50)^2]+exp[1.0],
    -sin[t]/t-1
};
Opt1[@f, optwaygrad,0, optnum,10];

结果（需多次运行）：

50.               50.               -1.151117991593894
45.04713670154168 50.73409357987649 -1.128374553525899
47.95358156651113 54.56960337517742 -1.128374553259392
48.61084727988142 45.18968474411513 -1.128374553040821
51.53888994360668 55.27956633578699 -1.11373603845343
50.0914003285785  38.652456006392   -1.070913459450462
55.57756844733133 40.1173971653779  -1.070913459450462
47.78264700546625 38.87080820770154 -1.070913459436176
43.59630987066741 40.6315049957481  -1.070913459402211
46.75204544703252 39.12675998582707 -1.070913459306036
10.

3 优化实例

3.1 较难的例子

3.1.1 拟合公式：y = (p1)+(p2*Exp(-p3*x/p5)+p4/(1+p4*p5*x))

p1,p2,p3,p4,p5为待求参数

数据（x， y）
0, 0.928
0.0000098, 1.02
0.0000195, 1.12
0.0000293, 1.25
0.0000391, 1.42
0.0000488, 1.7
0.0000586, 2.01
0.0000684, 2.26
0.0000781, 2.46
0.0000879, 2.63
0.0000977, 2.82
0.0001074, 3.01
0.0001172, 3.2
0.000127,  3.41
0.0001367, 3.59
0.0001465, 3.72
0.0001562, 3.85
0.000166,  3.98
0.0001758, 4.08

（正解的均方差RMSE＝0.033377163531，残差平方和为0.0211666658630，相关系数R ＝ 0.999497560）

Lu代码：

!!!using["luopt","math"];
f(p1,p2,p3,p4,p5:i,s,x,y:Array,max)=
{
    s=0,i=0,(i<max).while{
        x=Array[i,0], y=Array[i,1],
        s=s+[(p1)+(p2*exp(-p3*x/p5)+p4/(1+p4*p5*x))-y]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=19,
    Array=new[reals,max,2].SetArray{
        0, 0.928,
        0.0000098, 1.02,
        0.0000195, 1.12,
        0.0000293, 1.25,
        0.0000391, 1.42,
        0.0000488, 1.7,
        0.0000586, 2.01,
        0.0000684, 2.26,
        0.0000781, 2.46,
        0.0000879, 2.63,
        0.0000977, 2.82,
        0.0001074, 3.01,
        0.0001172, 3.2,
        0.000127,  3.41,
        0.0001367, 3.59,
        0.0001465, 3.72,
        0.0001562, 3.85,
        0.000166,  3.98,
        0.0001758, 4.08
    },
    Opt[@f, optdeep,50, optwaysimdeep, optwayconfra]
    //Opt1[@f, optdeep,5, optwaylm]    //也可以使用此语句
    //Opt1[@f, optwaynewton]  //也可以使用此语句
};

结果（需多次求解）：

6.855486033798995 4.813449005068769 -54298054.6899329 -10.72898361613976 -1516.502649185066 3.337716353170544e-002

3.1.2 拟合公式：z = p0*(1-exp(-p1*(x-p2)))+p3*x^p4+p5*x*y;

参数：p0 - p5
变量：x,y,z
数据(x,y,z)：
2        101       172
3        14        210
4        136       241
5        52        265
6        67        280
7        81        289
8        54        294
9        20        302
10        6        299
11        2        306

Lu代码：

!!!using["luopt","math"];
f(p0, p1, p2, p3, p4, p5 :i,s,x,y,z:Array,max)=
{
    s=0,i=0,(i<max).while{
        x=Array[i,0], y=Array[i,1], z=Array[i,2],
        s=s+[ p0*(1-exp(-p1*(x-p2)))+p3*x^p4+p5*x*y-z]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=10,
    Array=new[real_s,max,3].SetArray{
        "2 101 172
        3 14 210
        4 136 241
        5 52 265
        6 67 280
        7 81 289
        8 54 294
        9 20 302
        10 6 299
        11 2 306"
    },
    Opt1[@f, optmode,20, optwaysimdeep, optwayconfra]
    //Opt1[@f, optmode,20, optwaynewton]  //也可以使用此语句
};

结果（需多次求解）：

-317.2342157576544 -0.1348786321078753 -1.566716633721257 -4.949558747880146 2.255435954088224 -3.24860876032662e-003 1.399181338701023

3.348594116944429e+016 -4.19109191708201e-015 6.622514890426836e-003 257.3949699136492 0.8214692578887788 -6.209917030973698e-003 1.265994908321251

3.1.3 拟合公式：y=a*(x-d)^4+b*(x-d)^2+c;

x          y
-0.08    20.26008
-0.065   19.72613
-0.05    19.501619
-0.03    18.72662
-0.015   18.58769
0.015    18.592199
0.03     18.88372
0.05     19.5453
0.065    19.88743
0.08     20.9914
0        18.12336

Lu代码：

!!!using["luopt","math"];
f(a, b, c, d :i,s,x,y:Array,max)=
{
    s=0,i=0,(i<max).while{
        x=Array[i,0], y=Array[i,1],
        s=s+[ a*(x-d)^4+b*(x-d)^2+c-y]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=11,
    Array=new[reals,max,2].SetArray{
    "-0.08 20.26008
    -0.065 19.72613
    -0.05 19.501619
    -0.03 18.72662
    -0.015 18.58769
    0.015 18.592199
    0.03 18.88372
    0.05 19.5453
    0.065 19.88743
    0.08 20.9914
    0 18.12336"
    },
    Opt[@f, optwaysimdeep, optwayconfra]
    //Opt[@f, optwaylm]    //也可以使用此语句
    //Opt[@f, optwaynewton]  //也可以使用此语句
    //Opt1[@f, optdeep,10, optwaylm]    //也可以使用此语句
    //Opt1[@f, optwaynewton]  //也可以使用此语句
};

结果（需多次求解）：

420.0265769742538 -170.7483372687845 35.81768548688493 -0.449622975242104 0.1543213604571683

3.1.4 已知参数方程如下：

    x=r*(cos(t)+(t-b)*sin(t))+x0
    y=r*(sin(t)-(t-b)*cos(t))+y0

    拟合数据：

x 15.5910 24.6601 33.3732 49.3445 62.7992 68.3962 73.1584 79.9955 83.0531
y 86.2142 83.7114 80.2618 70.7216 58.0891 50.8058 42.9946 26.1718 8.4225

    中间变量 t 未知，待求参数 x0、y0、r 和 b。

Lu代码：

!!!using["luopt","math"];
f(x0, y0, r, b, t0, t1, t2, t3, t4, t5, t6, t7, t8 : i,s,t : tt,x,y,max)=
{
    tt.SetArray[t0, t1, t2, t3, t4, t5, t6, t7, t8],
    s=0,i=0,(i<max).while{
        t=tt[i],
        s = s + [ r*(cos(t)+(t-b)*sin(t))+x0 - x(i)]^2 + [ r*(sin(t)-(t-b)*cos(t))+y0 - y(i)]^2,
        i++
    },
    s
};
main(::tt,x,y,max)=
{
    max=9,
    tt=new[reals,max],
    x=new[reals,max].SetArray{"15.5910 24.6601 33.3732 49.3445 62.7992 68.3962 73.1584 79.9955 83.0531"},
    y=new[reals,max].SetArray{"86.2142 83.7114 80.2618 70.7216 58.0891 50.8058 42.9946 26.1718 8.4225"},
    Opt1[@f, optdeep,5, optwaysimdeep, optwayconfra]
};

结果（x0, y0, r, b, t0, t1, t2, t3, t4, t5, t6, t7, t8, 最小值）：

-3.296614781814681e-004 3.746723959397832e-004 -3.199522732021389 -27.5800929632305 -0.215438422526334 -0.3231575943836935 -0.4308780585413404 -0.646316754999914 -0.8617567278697795 -0.9694759174586836 -1.077195660800947 -1.292634766151454 -1.508072178584699 3.258061337614307e-009

3.2 非线性拟合测试题

参考：http://www.7d-soft.com/Regression_Test.htm

第1题：

!!!using["luopt","math"];
f(p1,p2,p3,p4,p5:i,s,x,y:Array,max)=
{
  s=0,i=0,(i<max).while{
    x=Array[i,0], y=Array[i,1],
    s=s+[1/(p1+p2*x^p3)+p4*x^p5-y]^2,
    i++
  },
  sqrt[s/max]
};
main(::Array,max)=
{
    max=178,
    Array=new[reals,max,2].SetArray{
160.73, 6266.7,
159.82, 6151.9,
158.84, 6035.1,
157.86, 5920.9,
156.87, 5812.6,
155.88, 5702.2,
154.89, 5594.9,
153.96, 5491.3,
152.97, 5385,
151.98, 5282.2,
150.99, 5181.3,
150.06, 5084.8,
149.08, 4988.8,
148.09, 4892.2,
147.1 , 4796.9,
146.17, 4701,
145.18, 4608,
144.2 , 4515.2,
143.2 , 4429.6,
142.21, 4342.9,
141.25, 4255.6,
140.2 , 4167.1,
139.14, 4077.6,
138.05, 3987.9,
136.96, 3898.9,
135.94, 3808.5,
134.84, 3717.7,
133.74, 3628.9,
132.65, 3543,
131.57, 3456.3,

130.55, 3372.5,
129.47, 3292.8,
128.4 , 3212.7,
127.33, 3133.6,
126.34, 3056.6,
125.29, 2985.5,
124.26 , 2912.5,
123.23, 2842.9,
122.21, 2774.1,
121.21, 2708.3,
120.27, 2642.1,
119.27, 2580.2,
118.29, 2518.7,
117.32, 2459.1,
116.42, 2401.1,
115.48, 2344.3,
114.55, 2290.9,
113.62, 2237.5,
112.7 , 2189,
111.85, 2138.8,
110.94, 2089.4,
110.03 , 2042.4,
109.13 , 1998.1,
108.28, 1953.6,
107.38, 1906.6,
106.48, 1867.8,
105.6 , 1824.5,
104.72 , 1784.3,
103.91, 1745,
103.05, 1704.5,

102.2, 1668.7,
101.35, 1629 ,
100.51, 1590.4,
99.739, 1552.5,
98.913 , 1514.7,
98.103, 1476.4,
97.308, 1444.3,
96.513, 1411.4,
95.78 , 1378.5,
95.002, 1344.8,
94.239, 1307.8,
93.482, 1276.1,
92.776, 1243.5,
92.039, 1212.9,
91.314, 1178.7,
90.604 , 1148.4,
89.942, 1115.9,
89.244 , 1084.5,
88.559, 1051.5,
87.889, 1029.7,
87.226, 996.16,
86.569 , 965.86,
85.963, 937.72,
85.323, 907.87,
84.694, 877.58,
84.081 , 838.17,
83.473, 819.48,
82.876, 797.76,
82.287, 768.54,
81.811, 749.96,

81.178, 724.39,
80.614 , 697.24,
80.118, 674.67,
79.574, 649.49,
79.011, 629.83,
78.478, 614.6,
78.012, 591.87,
77.494 , 573.43,
76.927, 558.94,
76.512, 539.45,
75.962, 526.99,
75.472, 514.14,
75.014, 504.11,
74.566, 484.4,
74.123, 473.23,
73.608, 468.93,
73.183, 453.77,
72.774, 448.58,
72.369, 447.73,
71.897, 431.79,
71.503, 432.45,
71.116, 432.15,
70.741, 420.71,
70.3, 427.26,
69.935, 419.76,
69.572, 407.28,
69.148, 408.04,
68.796, 393.71,
68.448, 403.74,
68.114, 408.8,

67.717, 401.26,
67.374, 400.81,
67.037, 401.89,
66.741, 408.68,
66.416, 398.49,
66.015, 414.14,
65.373, 419.78,
64.769, 426.82,
64.109, 418.42,
63.44 , 446.32,
62.772, 451.55,
62.111, 473.27,
61.508, 499.69,
60.908, 523.66,
60.219, 551.47,
59.699, 593.53,
59.119, 608.69,
58.547, 658.08,
57.992, 712.27,
57.483, 769.4,
56.969, 826.48,
56.472, 896.05,
55.989, 957.57,
55.513, 1065.1,
55.088, 1114.1,
54.651, 1195,
54.237, 1271.5,
53.836, 1355.6,
53.318, 1483.2,
52.701, 1690,

52.08 , 2245.9,
51.431, 2470.4,
50.877, 2719.1,
50.298, 2957.5,
49.74 , 3155.2,
49.2 , 3279.4,
48.702, 3546.4,
48.182, 3741,
47.681, 4021,
47.213, 4015.1,
46.768, 4304.7,
46.368, 4127.9,
45.956, 4530.9,
45.55 , 4802.9,
45.157, 5047.4,
44.799, 4804.3,
44.43 , 5164.1,
44.078, 4781,
43.727, 5175.5,
43.384, 5708.6,
43.079, 5679.6,
42.899, 5161.8,
42.719, 5399.1,
42.547, 5483,
42.253, 4839.4,
41.962, 5360.7,
41.691, 5622,
40.602, 5772.3
    },
    Opt[@f, optdeep,15]
    //Opt1[@f, optdeep,5]
};

结果（需尝试几次）：

1.773989909927714e-004 7.12634775546666e-033 16.72488443285935 1.215887381574098e-003 3.043657279780005 104.376667977067

第2题：

!!!using["luopt","math"];
f(p1,p2,p3,p4,p5,a1,a2,a3,a4:i,s,x1,x2,x3,x4,y:Array,max)=
{
    s=0,i=0,(i<max).while{
        x1=Array[i,0], x2=Array[i,1], x3=Array[i,2], x4=Array[i,3], y=Array[i,4],
        s=s+[(p1+p2*x1+p3*x2+p4*x3+p5*x4)/(1+a1*x1+a2*x2+a3*x3+a4*x4)-y]^2,
        i++
    },
    sqrt[s/max]
};
init(::Array,max)=
{
    max=18,
    Array=new[reals,max,5].SetArray{
15100, 29000, 508.0, 180, 3.40,
20500, 43350, 453.7, 141, 3.00,
80000, 92610, 487.9, 132, 2.70 ,
91500, 142775, 572.3, 182, 3.37,
82500, 2123160, 455.7, 113, 6.89 ,
20000, 227800, 481.3, 170, 5.03 ,
17800, 140000, 541.3, 179, 3.55 ,
3900,  15980, 538.6, 186, 2.72 ,
17300, 223200, 460.6, 100, 4.05 ,
25700, 229400, 393.1, 133, 3.22 ,
49400, 424500, 373.9, 106, 2.65 ,
40700, 561700, 482.8, 107, 1.91 ,
77000, 563600, 482.1, 140, 3.00 ,
92900, 557600, 415.1, 121, 1.31 ,
63300, 528300, 536.7, 144, 2.33 ,
51600, 488940, 385.1, 154, 3.55 ,
60000, 480500, 412.2, 111, 3.37 ,
70000, 530500, 567.2, 139, 2.55
    },
    Opt1[@f, optmax,1000, optmode,20, optdeep,20, optwaysimdeep, optwaylm, optwaygrad,0, optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5]
    //Opt[@f, optmax,1000, optmode,20, optdeep,20, optwaysimdeep, optwaylm, optwaygrad,0, optrange,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5,-1e5,1e5]
};

结果：

2.998714955530235 -6.9215674411842e-006 2.861867663835188e-007 -6.240030390485511e-003 4.853662624170808e-004 -2.235227336615253e-006 7.265550766777909e-008 -2.12622141557198e-003 3.084174242851936e-004 0.3028217094934776

第3题：

!!!using["luopt"];
g(x,y,p1,p2,p3)={p1/(1+p2/x+x/p3)-y}^2;
f(p1,p2,p3)=
{
    sqrt{[g(80.0 , 6.64 ,p1,p2,p3)+
    g(140.9 , 11.54 ,p1,p2,p3)+
    g(204.7 , 15.89 ,p1,p2,p3)+
    g(277.9 , 20.16 ,p1,p2,p3)+
    g(356.8 , 21.56 ,p1,p2,p3)+
    g(453.0 , 21.69 ,p1,p2,p3)+
    g(505.6 , 22.66 ,p1,p2,p3)+
    g(674.5 , 23.15 ,p1,p2,p3)+
    g(802.32 , 18.16 ,p1,p2,p3)+
    g(936.04 , 16.81 ,p1,p2,p3)]/10}
};
Opt[@f];
//Opt1[@f]  //也可以使用此语句

结果（需尝试几次）：

-101.0805795116041 -1258.521466815511 -170.1106261401753 0.8767278704475284

第4题：

!!!using["luopt","math"];
f(a0, a1, a2, a3, a4,b1,b2,b3,b4:i,s,x1,x2,x3,x4,y:Array,max)=
{
    s=0,i=0,(i<max).while{
        x1=Array[i,0], x2=Array[i,1], x3=Array[i,2], x4=Array[i,3], y=Array[i,4],
        s=s+[(a0+a1*x1+a2*x2+a3*x3+a4*x4)/(1+b1*x1+b2*x2+b3*x3+b4*x4)-y]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=23,
    Array=new[reals,max,5].SetArray{
        14 , 1.38 , -34, 16, 582 ,
        10, 0.52 , -29, 2 , 458 ,
        13 , 1.70 , -32 , 13 , 559 ,
        24 , 0.80 , 24 , 1 , 322 ,
        12 , 1.83 , 41 , 11, 399 ,
        6 , 1.77 , -50 , 7 , 523 ,
        18, 1.23 , 27 , 4 , 322 ,
        -10, 0.28, -8, 6 , 358 ,
        0 , 1.20 , 66 , 6 , 354 ,
        14 , 1.75, -60 , 6 , 574 ,
        12 , 1.78 , -70 , 7 , 489 ,
        -18, 1.37, -15 , 0 , 232 ,
        16, 1.38, 0 , 4 , 440 ,
        -4, 0.29, -9 , -7, 421 ,
        -23, 1.12, -12, -14, 181 ,
        5 , 1.52, 0 , 10 , 426 ,
        -16, 0.63, 34 , 1 , 364 ,
        -1 , 1.32 , 22 , -7 , 375 ,
        -18, 1.18, 4 , -11, 224 ,
        8., 1.50 , -11 , 5 , 514 ,
        -8 , 1.43, 4 , -12, 381 ,
        -11, 0.74, 10 , 0 , 275 ,
        -19, 1.07, -5, 0 , 426
    },
    Opt1[@f, optmax,1000, optmode,20, optdeep,20, optwaysimdeep, optwayconfra, optrange : -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10]
    //Opt1[@f, optmax,1000, optmode,20, optwaygrow, optrange : -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10]  //也可以使用此语句
    //Opt[@f, optmax,1000, optmode,20, optdeep,20, optwaysimdeep, optwayconfra, optrange : -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10]  //也可以使用此语句
};

结果（需求解几次）：

674.6781937042127 227.7453915656436 2120.317418232974 1.643200016850861 -176.0515628702959 0.5725832815136416 5.556417489954765 3.344065962703589e-002 -0.5600170829989317 48.05714060993065

第5题：

!!!using["luopt","math"];
f(p1, p2, p3, p4, p5,p6,p7,p8:i,s,x,y,z:Array,max)=
{
    s=0,i=0,(i<max).while{
        x=Array[i,0], y=Array[i,1], z=Array[i,2],
        s=s+[ p1+p2*x^p3+p4*y^p5+p6*x^p7*y^p8-z]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=24,
    Array=new[reals,max,3].SetArray{
500, 25, 1.5 ,
500, 50, 2.25 ,
500, 100 ,3.15 ,
500, 200, 4.0 ,
500, 300, 4.2 ,
500, 400, 4.3 ,
1000, 25 , 1.45 ,
1000, 50 , 2.35 ,
1000, 100, 3.95 ,
1000, 200, 6.95 ,
1000, 300, 8.15 ,
1000, 400, 8.4 ,
1500, 25 , 1.45 ,
1500, 50 , 2.45 ,
1500, 100, 4.15 ,
1500, 200, 7.45 ,
1500, 300, 10.65 ,
1500, 400, 11.85 ,
2000, 25 , 1.45 ,
2000, 50 , 2.5 ,
2000, 100, 4.2 ,
2000, 200, 7.75 ,
2000, 300, 11.45 ,
2000, 400, 14.3
    },
    Opt1[@f, optmax,500, optmode,20, optdeep,20, optwaysimdeep, optwayconfra]
    //Opt[@f, optmax,500, optmode,50, optdeep,50, optwaysimdeep, optwayconfra]  //也可以使用此语句
};

结果（需多次求解）：

第6题：

!!!using["luopt","math"];
f(a0, a1, a2, a3, k1,k2,k3 :i,s,x,y:Array,max)=
{
    s=0,i=0,(i<max).while{
        x=Array[i,0], y=Array[i,1],
        s=s+[ a0+a1*x^k1+a2*x^k2+a3*x^k3-y]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=10,
    Array=new[reals,max,2].SetArray{
        1.0, 8.2,
        2.0, 4.6,
        3.0, 4.3,
        4.0, 4.6,
        5.0, 5.1,
        6.0, 5.5,
        7.0, 5.7,
        8.0, 5.5,
        9.0, 5.0,
        10.0,3.8
    },
    Opt1[@f, optwaysimdeep, optwayconfra]
    //Opt[@f, optwaysimdeep, optwayconfra]  //也可以使用此语句
};

结果（需多次求解）：

-2.435047968667117 -5.596012968756312e-004 1.746941301238457 8.888535605114285 4.021264991475747 0.8188562111431067 -1.1645140483289 2.147261368496635e-002

第7题：

!!!using["luopt","math"];
f(p1, p2, p3, p4, p5,p6,p7:i,s,x,y,z:Array,max)=
{
    s=0,i=0,(i<max).while{
        x=Array[i,0], y=Array[i,1], z=Array[i,2],
        s=s+[ (p1+p2*x+p3*y+p4*x*y)/(1+p5*x+p6*y+p7*x*y)-z]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=21,
    Array=new[reals,max,3].SetArray{
4332, 26.94 , 43.70 ,
4697, 23.64 , 44.50 ,
5062, 25.19 , 47.70 ,
5428, 28.60 , 52.30 ,
5793, 28.74 , 54.21 ,
6158, 29.33 , 55.58 ,
6523, 32.92 , 55.70 ,
6889, 31.87 , 57.70 ,
7254, 33.07 , 58.60 ,
7619 , 29.36 , 60.24 ,
7984, 27.96 , 59.13 ,
8350 , 30.18 , 57.00 ,
8715 , 37.84 , 57.30 ,
9080 , 38.86 , 59.00 ,
9445 , 35.18 , 60.20 ,
9811 , 28.81 , 61.80 ,
10176, 34.57 , 63.17 ,
10541, 37.49 , 66.23 ,
10906, 29.30 , 61.80 ,
11272, 32.47 , 62.03 ,
11637, 38.24 , 65.30
    },
    Opt1[@f, optmax,1000, optmode,20, optdeep,20, optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e5,1e5,-1e5,1e5]
    //Opt[@f, optmax,1000, optmode,200, optdeep,50, optwaysimdeep, optwayconfra, optrange,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e10,1e10,-1e5,1e5,-1e5,1e5]
};

结果（需多次求解）：

92.07386083215749 -2.673470802502084e-002 -2.720784311480835 7.444462298472703e-004 -3.845503645661088e-004 -3.039200875467085e-002 1.070399375206463e-005 1.006260685262059

第8题：

!!!using["luopt","math"];
f(p1, p2, p3, p4, p5,p6:i,s,x1,x2,x3,y:Array,max)= //目标函数定义
{
    s=0,i=0,(i<max).while{
        x1=Array[i,0], x2=Array[i,1], x3=Array[i,2], y=Array[i,3],
        s=s+[ p1/((p2+x1)*(1+p3*x2)*(x3-p4)^2)+p5*x3^p6-y]^2,
        i++
    },
    sqrt[s/max]
};
main(::Array,max)=
{
    max=27,
    Array=new{reals,max,4, data:
        34.9, 0.043378, 8., 0.996556 ,
        34.9, 0.216888, 8., 0.985708 ,
        34.9, 0.433776, 8., 0.973826 ,
        58.2, 0.026027, 8., 0.999409 ,
        58.2, 0.130133, 8., 0.99817 ,
        58.2, 0.260265, 8., 1.000176 ,
        93.1, 0.016267, 8., 0.995131 ,
        93.1, 0.081333, 8., 1.009887 ,
        93.1, 0.162666, 8., 1.008251 ,
        34.9, 0.043378, 20., 0.835576 ,
        34.9, 0.216888, 20., 0.777734 ,
        34.9, 0.433776, 20., 0.715483 ,
        58.2, 0.026027, 20., 0.854949 ,
        58.2, 0.130133, 20., 0.822743 ,
        58.2, 0.260265, 20., 0.784273 ,
        93.1, 0.016267, 20., 0.85902 ,
        93.1, 0.081333, 20., 0.841512 ,
        93.1, 0.162666, 20., 0.81895 ,
        34.9, 0.043378, 40., 0.387322 ,
        34.9, 0.216888, 40., 0.338941 ,
        34.9, 0.433776, 40., 0.293558 ,
        58.2, 0.026027, 40., 0.342388 ,
        58.2, 0.130133, 40., 0.311761 ,
        58.2, 0.260265, 40., 0.280112 ,
        93.1, 0.016267, 40., 0.308071 ,
        93.1, 0.081333 , 40., 0.287257 ,
        93.1, 0.162666, 40., 0.264443
    },
    Opt1[@f, optmode,9, optdeep,9, optwaysimdeep, optwayconfra]
    //Opt[@f, optwaysimdeep, optwayconfra]  //也可以使用此语句
};

结果（需多次求解）：

178963.6676067503 3672.74944952677 0.530211167053818 27.80500822839966 195.3596009106615 -2.596554487779508 1.977698395027922e-002

第9题：

!!!using["luopt"];
y(x,a,b,c,d)=a*exp(b*abs(x+c)^d);
f(a,b,c,d)=
{
    sqrt{{[y(27.25 ,a,b,c,d)-1]^2
    +[y(27.75 ,a,b,c,d)-3]^2
    +[y(28.25 ,a,b,c,d)-6]^2
    +[y(28.75 ,a,b,c,d)-13]^2
    +[y(29.25 ,a,b,c,d)-18]^2
    +[y(29.75 ,a,b,c,d)-19]^2
    +[y(30.25 ,a,b,c,d)-17]^2
    +[y(30.75 ,a,b,c,d)-16]^2
    +[y(31.25 ,a,b,c,d)-6]^2
    +[y(31.75 ,a,b,c,d)-5]^2
    +[y(32.25 ,a,b,c,d)-2]^2}/11}
};
Opt1[@f, optwaylm, optwaygrad,0];
//Opt1[@f, optwaysimdeep, optwayconfra];   //也可以使用此语句
//Opt[@f, optwaysimdeep, optwayconfra];   //也可以使用此语句
//Opt[@f, optwaygrow, optwaygrad,0];    //也可以使用此语句
　

结果（需多次求解）：

19.15817773649649 -0.3625927565554836 -29.81592272236692 2.297951055806539 1.15469090018263

3.3 约束优化

3.3.1求最小值 f(x,y)= -(x*sin(9*pi*y)+y*cos(25*pi*x)+20)，要求 x^2+y^2<=81。

Lu代码：

!!!using["luopt","math"];
f(x,y)= -(x*sin(9*pi*y)+y*cos(25*pi*x)+20);
g(x,y)= which[x^2+y^2<=81, 1,0];
Opt1[@f, optconfun, @g, optmax,200, optmode,10, optwaygrow, optwayconfra, optwaygrad,0];

结果：

-6.440025821695283 -6.277972013333297 -32.71788780688353

3.3.2 求最大值 (sin(2*pi*x1))^3*sin(2*pi*x2)/(x1^3*(x1+x2));

要求：x1^2-x2+1 <= 0;
    1-x1+(x2-4)^2 <= 0;
    0<=x1<=10;
    0<=x2<=10;

Lu代码：

!!!using["luopt","math"];
f(x1,x2)= - [(sin(2*pi*x1))^3*sin(2*pi*x2)/(x1^3*(x1+x2))];  //转换为求最小值
g(x1,x2)= which[x1^2-x2+1 <= 0 & 1-x1+(x2-4)^2 <= 0 : 1, 0];
Opt1[@f, optconfun, @g, optmax,200, optmode,100, optrange : 0.0,10.0; 0.0,10.0];

结果：

1.227970948659768 4.245370242657014 -9.582504139865343e-002

最大值为：9.582504139865343e-002

3.3.3 求最小值 10*x1+9*x2+8*x3+7*x4*sin(x1+x2+x3);

要求：(3*x2+2*x4*cos(x1+x2+x3+x4))^2<=90;
    x1+x2>=-30;
    x3+x4>=30;
    3*x1+2*x3<=120;

参数范围均在[-100,100]之间。

Lu代码：

!!!using["luopt","math"];
f(x1,x2,x3,x4)= 10*x1+9*x2+8*x3+7*x4*sin(x1+x2+x3);
g(x1,x2,x3,x4)= which[(3*x2+2*x4*cos(x1+x2+x3+x4))^2<=90 & x1+x2>=(-30) & x3+x4>=30 & 3*x1+2*x3<=120, 1,0];
Opt1[@f, optconfun, @g, optmax,200, optmode,200, optdeep,20, optrange : -100.0,100.0; -100.0,100.0; -100.0,100.0; -100.0,100.0];

结果：

-98.78718784979807 68.78718939245172 -65.84111673048275 99.09825760408047 -1589.027715572939

3.3.4 约束拟合 y = (p3+p2*p1*x^p4)/(p5+p1*x^p4);

要求：p2=3*p3+p1+p2+p3;
    20.3>=p1+p2>=20;

数据：

x：0.010, 0.020, 0.040, 0.060, 0.080, 0.100, 0.120, 0.140, 0.160, 0.180, 0.200, 0.220, 0.240, 0.260, 0.280, 0.300, 0.320, 0.340, 0.360, 0.380, 0.400
y：3.936, 4.117, 4.775, 5.553, 6.268, 6.935, 7.480, 8.188, 8.361, 8.533, 8.771, 9.120, 9.024, 9.288, 9.462, 9.379, 9.685, 9.482, 9.545, 9.604, 9.546

分析：本例中等式约束可以化简，本例按不化简拟合。等式约束可取平方并乘以一个较大的数加到目标函数中。

Lu代码：

!!!using["luopt","math"];
g(p1,p2,p3,p4,p5)= which[20.3>=p1+p2 & p1+p2>=20 : 1,0];
f(p1,p2,p3,p4,p5 : i,s,x,y : xx,yy,max)=
{
    i=-1,s=0,while{++i<max,
    x=xx[i], y=yy[i],
    s=s+[(p3+p2*p1*x^p4)/(p5+p1*x^p4) - y]^2 + 1e10*[3*p3+p1+p2+p3 - p2]^2 // 注意等式约束添加到目标函数的方法
    },
    s
};
main(:: xx,yy,max)=
{
    xx=ra1[0.010, 0.020, 0.040, 0.060, 0.080, 0.100, 0.120, 0.140, 0.160, 0.180, 0.200, 0.220, 0.240, 0.260, 0.280, 0.300, 0.320, 0.340, 0.360, 0.380, 0.400],
    yy=ra1[3.936, 4.117, 4.775, 5.553, 6.268, 6.935, 7.480, 8.188, 8.361, 8.533, 8.771, 9.120, 9.024, 9.288, 9.462, 9.379, 9.685, 9.482, 9.545, 9.604, 9.546],
    max=len[xx],
    Opt1[@f, optconfun, @g, optmax,500, optmode,200, optdeep,20]
};

结果：

1.174560633533379e-014 20.0000000733822 7.391404517515066e-015 0.4861658465496704 8.143677889916323e-015 2.235330571817288

3.4 多维优化

3.4.1 求下列多维函数最小值： ∑(3*(cos(2*x[i]) + sin(2*x[i+1])) + sqrt(x[i+1]^2 + x[i]^2));

    其中：n=20，i=0～n-2，x[i]∈[-30,30]。

Lu代码：

!!!using["luopt","math"];
f(a : i, s , max)=
{
    max=len[a]-1,
    s=0, i=0, (i<max).while{
        s=s+(3*(cos(2*a[i]) + sin(2*a[i+1])) + sqrt(a[i+1]^2 + a[i]^2)),
        i++
    },
    s
};
main(:Array, max)=
{
    max=20,
    Array=new[real_s, max],
    Opt1[Array, @f, optrange : -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0; -30.0,30.0]
};

结果：-51.7695459159333

1.503381585935963 -1.101177339173767 -1.094504428145357 -1.094373225161253 -1.094328127754031 -1.094326554999183 -1.094445477955365 -1.094357292227312 -1.094410614814496 -1.094385668267113 -1.094393836108593 -1.094363347681921 -1.094375183173841 -1.094251791948202 -1.09434566811177 -1.094401404833669 -1.09436516172288 -1.094313702765746 -1.087307395163289 -0.7385347279659271 -51.7695459159333

3.4.2 求下列多维函数最小值： ∑(3*(cos(2*x[i]) + sin(2*x[i+1])) + sqrt(x[i+1]^2 + x[i]^2));

    其中：n=100，i=0～n-2。

Lu代码：

!!!using["luopt","math"];
f(a : i, s , max)=
{
    max=len[a]-1,
    s=0, i=0, (i<max).while{
        s=s+(3*(cos(2*a[i]) + sin(2*a[i+1])) + sqrt(a[i+1]^2 + a[i]^2)),
        i++
    },
    s
};
Opt1[new(real_s, 100), @f, optmax,200, optmode,10, optdeep,10];  //n=100，无约束优化

结果：-262.6193574331547

-1.503589794995611 -1.101165993790981 -1.094485973512971 -1.094477365687782 -1.094348063147205 -1.094374887930401 -1.094406810150064 -1.094254437325214 -1.094354679300382 -1.094357945608955 -1.094386693983307 -1.094400195363591 -1.094430881361193 -1.094301980400329 -1.094451983562294 -1.09440057303286 -1.094246391707527 -1.094393611497804 -1.094317884914888 -1.094331829529994 -1.094416479129036 -1.094398963357002 -1.094409209599468 -1.094269021137714 -1.09440364315413 -1.09448075055456 -1.094345890346567 -1.094377882486256 -1.094335811282484 -1.094311557019327 -1.094181460610167 -1.094355995795314 -1.094362619422068 -1.094391693363947 -1.094396242056888 -1.094389706369454 -1.094308176206261 -1.094236009983625 -1.094468152806394 -1.094264582930009 -1.094218568945094 -1.094487868407586 -1.094371004727736 -1.094335754636069 -1.094493710249255 -1.094413947834062 -1.094395743895835 -1.094460063524959 -1.094315544253196 -1.094309224794404 -1.094238853251452 -1.094299581429759 -1.094460738638513 -1.094440093985098 -1.094384292459407 -1.094465833693187 -1.094300044371981 -1.094333544439325 -1.094372253887107 -1.094371552346044 -1.094412279004186 -1.094350558799062 -1.094320418874927 -1.094364991752822 -1.094641012221824 -1.094396271927086 -1.094332090889283 -1.094255353249292 -1.094441019469457 -1.094386255277219 -1.094304210929158 -1.0942951145132 -1.09439223657837 -1.094371933599244 -1.09444717731441 -1.094369073964514 -1.094269805994881 -1.094449953762979 -1.094352352398937 -1.09442676376325 -1.094302907367715 -1.094415111889624 -1.094237968076476 -1.094374925733766 -1.094444387927041 -1.094394454068458 -1.094383346719542 -1.094219317247797 -1.094300421296237 -1.09434971557079 -1.094417944826193 -1.094373380532534 -1.09433164816419 -1.094268795813747 -1.094479931551439 -1.094322348615986 -1.094382411390265 -1.094152442923652 -1.087374832951749 -0.7384077798609477 -262.6193574331547

3.4.3 求下列多维函数（Extended Freudenstein & Roth Function）最小值： ∑[-13+x(2*i-2)+((5-x(2*i-1))*x(2*i-1)-2)*x(2*i-1)]^2 + [-29+x(2*i-2)+((x(2*i-1)+1)*x(2*i-1)-14)*x(2*i-1)]^2;

    其中：n=2,4,6,8,... ...（n取2的倍数），i=1～n/2，拟合参数x[j]中j=0～n-1，共n个拟合参数。理论最小值为0，最优参数为 5,4,5,4, ... ...。

分析：本例随n的增大，拟合难度显著增大，耗时较长。

Lu代码：

!!!using["luopt","math"];
f(x : i ,s,n) =
{
    n=len(x)/2, i=0, s=0, while{++i<=n, s=s+[-13+x(2*i-2)+((5-x(2*i-1))*x(2*i-1)-2)*x(2*i-1)]^2 + [-29+x(2*i-2)+((x(2*i-1)+1)*x(2*i-1)-14)*x(2*i-1)]^2 } , s
};
main(:Array, n)=
{
    n=20,  //n从2开始，分别取2的倍数，观察最小值输出窗口的中间结果输出情况。可以取更大的n试试，例如 n = 50 甚至更大的数值。
    Array=new[real_s, n],
    Opt1[Array, @f, optmax,1000, optmode,20, optdeep,21]
    //Opt1[Array, @f, optmax,1000, optmode,20, optdeep,21, optwaycom, optwaylm, optwaynewton, optwaysimdeep, optwayconfra, optwaygrad,0]  //也可以使用此语句
};

结果：

4.999999999999984 4. 4.999999999999998 4. 5.000000000000005 4. 4.999999999999993 4. 5.000000000000002 4. 5.000000000000006 4. 4.999999999999993 4. 4.999999999999998 4. 5.000000000000004 4. 5.000000000000002 4. 8.6538041302745e-028

3.4.4 求下列多维函数（Extended Miele & Cantrell Function）最小值： ∑exp[x(4*i-4)-x(4*i-3)]^2 + 100*[x(4*i-3)-x(4*i-2)]^6 + [tan(x(4*i-2)-x(4*i-1))]^4 + x(4*i-4)^8;

    其中：n=4,8,16,20,... ...（n取4的倍数），i=1～n/4，拟合参数x[j]中j=0～n-1，共n个拟合参数。理论最小值为0，最优参数为 0,1,1,1, 0,1,1,1, ... ...。

分析：本例仅作演示用。由于计算精度所限（使用双精度实数，16位有效数字），本例的最优参数难以得到，但可以使最小值接近于0。

Lu代码：

!!!using["luopt","math"];
f(x : i ,s,n) =
{
    n=len(x)/4, i=0, s=0, while{++i<=n, s=s+exp[x(4*i-4)-x(4*i-3)]^2 + 100*[x(4*i-3)-x(4*i-2)]^6 + [tan(x(4*i-2)-x(4*i-1))]^4 + x(4*i-4)^8} , s
};
main(:Array, n)=
{
    n=20,  //n从4开始，分别取4的倍数，观察最小值输出窗口的中间结果输出情况。可以取更大的n试试，例如 n=100。
    Array=new[real_s, n],
    Opt1[Array, @f, optmax,200, optwaycom, optwaysimdeep, optwayconfra]
};

4 求方程（组）的全部解

[返回页首] luopt::iFind(f, luopt::optmax,m, luopt::optrange,min,max, luopt::opteps,eps, luopt::optexact,exact, luopt::optpara,x1,x2,..., luopt::optthis,i, luopt::optmin0,min0, luopt::optsameeps,sameeps, luopt::optsameabs,sameabs, luopt::optstarter,starter, luopt::optarray,a, luopt::optarraybegin,&k, luopt::optout,out, luopt::optnum,num)：求单变量方程的全部解

f：自定义n元函数句柄，不可缺省。格式如下：

f(x1,x2,...,xn)=
{
    ... ...
};

    默认求方程f(x1,x2,...,xn)=0第一个自变量的全部解，而其他自变量赋值为0.0。可以用参数optthis指定求解的自变量，也可以用参数optpara给出其他自变量的值。

luopt::optmax,m：区间分割数目（大于等于10），区间分割数目越多精度越高。可以缺省该参数，缺省值为200。
luopt::optrange,min,max：指定求解区间。若缺省该参数，则min为-1e50，max为1e50。
luopt::opteps,eps：控制精度要求，eps>0。可以缺省该参数，缺省值为1e-6。
luopt::optexact,exact：控制精度的方式，exact非0将严格控制精度，否则只追求解的正确性，尽可能满足解的精度。可以缺省该参数，将严格控制解的精度。
luopt::optpara,x1,x2,...：给指定求解自变量之外的其他自变量赋值，参数x1,x2,...的个数比全部自变量个数少1。若缺省该参数，其他自变量缺省值均为0.0。
luopt::optthis,i：指定求解的自变量。0表示第一个自变量；1表示第二个自变量，以此类推。若缺省该参数，i为0。
luopt::optmin0,min0：解的可能的极小值。可以缺省该参数，缺省值为±1e-50。
luopt::optsameeps,sameeps：相对比较两组值是否相等的极小值。sameeps>0。可以缺省该参数，缺省值为1e-3。
luopt::optsameabs,sameabs：绝对比较两组值是否相等的极小值。sameeps>=0。可以缺省该参数，缺省值为0.0。
luopt::optstarter,starter：提供一个初值。可以缺省该参数。
luopt::optarray,a：提供一个数组，用于接受结果。一组结果包括所有的自变量参数及误差。默认从数组开始存放数据，也可用参数optarraybegin指定开始存放数据的位置。可以缺省该参数。
luopt::optarraybegin,&k：指定数组中开始存放数据的位置，k返回下一个存放位置。可以缺省该参数，默认从地址0开始存放。
luopt::optout,out：是否输出结果。out非0时输出。可以缺省该参数，默认输出结果。
luopt::optnum,num：设置最多求解的解的个数。可以缺省该参数，默认最多求解100个解。

返回值：解的个数。

运行错误：

1：指定的表达式不存在；
2：常量表达式；
3：自变量参数重新赋值，这是不允许的；
4：参数不匹配；
5：不可识别描述符；
6：参数不符合要求；
7：自定义函数返回了非法值。

例子1：求方程的全部解：f(x)=x^6-5*x^5+3*x^4+x^3-7*x^2+7*x-20;

f(x)=x^6-5*x^5+3*x^4+x^3-7*x^2+7*x-20;
luopt::iFind[@f];

例子2：求方程的全部解（x取-8～8，y取1）：f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;

!!!using("luopt");
f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;
iFind[@f,optrange,-8.,8.,optpara,1.];

例子3：求方程的全部解（x取1，y取-8～8）：f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;

!!!using("luopt");
f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;
iFind[@f, optrange,-8.,8., optthis,1, optpara,1.];

[返回页首] luopt::Find(f, luopt::optmax,m, luopt::optmode,mode, luopt::optdeep,deep, luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax, luopt::optautorange,autorange，luopt::opteps,eps, luopt::optmin0,min0, luopt::optsameeps,sameeps, luopt::optsameabs,sameabs, luopt::optstarter,&x1,&x2,...,&xn,&min, luopt::optarray,a, luopt::optarraybegin,&k, luopt::optexpand,expand, luopt::optcontract,contract, luopt::optdx,dx, luopt::optstep,step, luopt::optout,out, luopt::optnum,num)：求方程组的全部解

f：自定义2*n元函数句柄，不可缺省，必须由二级函数HFor获得该句柄。格式如下：

f(x1,x2,...,xn,y1,y2,...,yn)=
{
    y1= ...,
    y2= ...,
    ... ...,
    yn= ...
};

luopt::optmax,m：随机点的数目（大于等于10），点数越多精度越高。可以缺省该参数，缺省值为200。
luopt::optmode,mode：工作模式，取0、1、2、3、... ...。通常，工作模式取值越大，搜到的解越多，但耗时越长。若缺省该参数，工作模式取0。
luopt::optdeep,deep：搜索深度，取0、1、2、3、... ...。通常，搜索深度取值越大，搜到的解越多，但耗时越长。若缺省该参数，搜索深度取0。
luopt::optrange,x1min,x1max,x2min,x2max,...,xnmin,xnmax：指定求解区间。若缺省该参数，则所有变量区间均为[-1e50～1e50]。
luopt::optautorange,autorange：autorange!=0时将在指定搜索区间内自动调节合适的搜索区间，否则不进行调节。默认是不会自动调节搜索区间的。
luopt::opteps,eps：控制精度要求，eps>0。可以缺省该参数，缺省值为1e-6。满足sqrt[(y1*y1+y2*y2+...+yn*yn)/n]<eps。
luopt::optmin0,min0：解的可能的极小值。可以缺省该参数，缺省值为±1e-50。
luopt::optsameeps,sameeps：相对比较两组值是否相等的极小值。sameeps>0。可以缺省该参数，缺省值为1e-2。
luopt::optsameabs,sameabs：绝对比较两组值是否相等的极小值。sameeps>=0。可以缺省该参数，缺省值为1e-3。
luopt::optstarter,&x1,&x2,...,&xn,&min：前n个参数提供一组初值，并返回一组解；最后一个参数min返回该组解的误差。可以缺省该参数。
luopt::optarray,a：提供一个数组，用于接受结果。一组结果包括所有的自变量参数及误差。默认从数组开始存放数据，也可用参数optarraybegin指定开始存放数据的位置。可以缺省该参数。
luopt::optarraybegin,&k：指定数组中开始存放数据的位置，k返回下一个存放位置。可以缺省该参数，默认从地址0开始存放。
luopt::optexpand,expand：扩张系数。一般取1.2～2.0。可以缺省该参数，缺省值为1.6。
luopt::optcontract,contract：收缩系数。一般取0～1。可以缺省该参数，缺省值为0.5。
luopt::optdx,dx：搜索增量初值，dx>0。可以缺省该参数，缺省值为0.1。
luopt::optstep,step：搜索步长修正系数。step取0～1之间的数。可以缺省该参数，缺省值为0.1。
luopt::optout,out：是否输出结果。out非0时输出。可以缺省该参数，默认输出结果。
luopt::optnum,num：设置最多求解的解的个数。可以缺省该参数，默认最多求解10个解。

返回值：解的个数。

说明：需要多次运行该函数以获得需要的解。

运行错误：

1：指定的表达式不存在；
2：常量表达式或者参数不成对；
3：内存错误；
4：参数不匹配；
5：不可识别描述符；
6：参数不符合要求；
7：自定义函数返回了非法值。

例子1：解方程组：

(x-y)^2-3*(x-y) = 10
x^2+2*x*y+y^2 = 9

代码：

f(x,y,y1,y2)=
{
  y1=(x-y)^2-3*(x-y)-10,
  y2=x^2+2*x*y+y^2-9
};
luopt::Find[@f];

例子2：解方程组：

2*x1-x2^2-exp(-x1) = 0
-(x1^3)+x1*x2-exp(-x2) = 0

代码：

f(x1,x2,y1,y2)=
{
  y1=2*x1-x2^2-exp(-x1),
  y2=-(x1^3)+x1*x2-exp(-x2)
};
luopt::Find[@f];

例子3：解方程组：t取-7～7

-b*sin(a+6*t)+n-40.4945=0
-b*sin(a+7*t)+n-40.5696=0
-b*sin(a+8*t)+n-41.0443=0
-b*sin(a+9*t)+n-41.4190=0

代码：

!!!using["luopt"];
f(a,b,n,t,y1,y2,y3,y4)=
{
  y1=-b*sin(a+6*t)+n-40.4945,
  y2=-b*sin(a+7*t)+n-40.5696,
  y3=-b*sin(a+8*t)+n-41.0443,
  y4=-b*sin(a+9*t)+n-41.4190
};
Find[@f, optmode,5, optrange,-1e50,1e50,-1e50,1e50,-1e50,1e50,-7.,7.];

5 聚类

[返回页首] luopt::PeakCluster(data, type, PSErr, PS, PPErr, PP, max, min, it)：利用数据曲线尖峰形状和位置聚类，需调整参数获得最佳聚类结果

data：实数数据矩阵（m×n），m行n列，每一行为一条数据曲线。
type：整数数组（m×3），返回聚类结果，第一列为曲线类别（用自1开始的整数标识），第二列为曲线形状相似度百分数，第三列为尖峰位置准确度百分数。
PSErr：整数，曲线形状沿水平轴的允许误差（偏移量），PSErr>=0，PSErr越大允许误差越大。
PS：实数，取值范围为[0.0～1.0]，曲线形状最低相似度。
PPErr：整数，尖峰位置沿水平轴的允许误差（偏移量），PPErr>=0，PPErr越大允许误差越大。
PP：实数，取值范围为[0.0～1.0]，尖峰位置最低准确度。
max,min：整数，用于比较的极大值和极小值数目。max和min不可同时为0；若max<0表示取所有的极大值；若min<0表示取所有的极小值。
it：最大迭代次数。

返回值：整数。-2表示所有行中均无极值；-1表示其他各类错误；返回值大于0表示聚类成功，返回实际迭代次数。

运行错误：

1：参数不匹配；
2：参数不符合要求；
3：内存错误。

例子1：共有8条数据曲线，每条曲线55个数据

代码：

!!!using["luopt","math","win","sys"];
main(: Array, row, column, ia, i, f, s)=
{
    Array=matrix[GetFileName().readtxt()], //读取文件中的数据并转换为一个矩阵，一条曲线数据为一行
    len[Array,0,&row, &column], //用len函数取矩阵的行数和列数
    ia=new[ints, row, 3],     //存放聚类结果
    PeakCluster[Array,ia,0,0.6,0,0.5,-1,-1, 100], //尖峰法聚类
    f=fopen["D:\\a.txt","wb+"], //聚类结果写入文件a.txt
    fwrite[f,UTEXT,2], s=new[string,30],
    i=-1, (++i<row).while{
        s[0]=0, sprintf(s,"%d\t%d\t%d", ia[i,0], ia[i,1], ia[i,2]), fwrite[f,s], fwrite[f,"\r\n"]
    },
    fclose[f]
};

聚类结果：

例子2：2021 年高教社杯全国大学生数学建模竞赛 E题 中药材的鉴别 问题1

6 函数泰勒展开

[返回页首] luopt::FunTaylor(mode, m, n, a, d)：将任意函数按给定系数根据泰勒公式展开并计算

mode：整数，工作模式。mode=0返回泰勒展开式的系数个数；mode=1返回字符串形式的泰勒展开式；mode=2计算函数值。
m：整数，任意函数的参数个数（m>1）。
n：整数，泰勒公式级数（n>1）。
a：一维实数数组，存放泰勒展开式的系数。
d：实数数据矩阵（k×（m+1）），k行m+1列，每一行前m个数存放函数参数，最后一个数返回函数值。

返回值：整数。mode=0或mode=1返回泰勒展开式的系数个数；mode=2时正常返回1，出错时返回0。

运行错误：

1：参数不匹配；
2：参数不符合要求；
3：内存错误。

例子1：

!!!using["luopt","math"];; //使用命名空间
main(:a,b,k,i)=
{
    k=FunTaylor[0,3,4], //返回函数有3个参数时，泰勒公式4级展开的系数个数
    a=new[real_s,k],   //申请数组，存放泰勒展开式系数
    i=-1, (++i<k).while{a[i]=i+1.0}, //数组初始化
    FunTaylor[1,3,4,a], //返回字符串形式的泰勒展开式
    b=matrix{       //存放函数参数，最后一列用于存放函数值
    "
    0.5 1 2 0
    0.5 2 1 0
    0.5 3 0 0
    "
    },
    FunTaylor[2,3,4,a,b], //根据给定系数及函数参数计算函数值
    o[b]          //输出计算结果
};

结果：可以看出，有38个系数

x0=y0-(1.), x1=y1-(2.), x2=y2-(3.),
y=(4.)+(5.)*x0+(6.)*x1+(7.)*x2
  +(8.)*x0*x0+(9.)*x0*x1+(10.)*x0*x2+(11.)*x1*x1+(12.)*x1*x2+(13.)*x2*x2
  +(14.)*x0*x0*x0+(15.)*x0*x0*x1+(16.)*x0*x0*x2+(17.)*x0*x1*x1+(18.)*x0*x1*x2+(19.)*x0*x2*x2+(20.)*x1*x1*x1+(21.)*x1*x1*x2+(22.)*x1*x2*x2+(23.)*x2*x2*x2
  +(24.)*x0*x0*x0*x0+(25.)*x0*x0*x0*x1+(26.)*x0*x0*x0*x2+(27.)*x0*x0*x1*x1+(28.)*x0*x0*x1*x2+(29.)*x0*x0*x2*x2+(30.)*x0*x1*x1*x1+(31.)*x0*x1*x1*x2+(32.)*x0*x1*x2*x2+(33.)*x0*x2*x2*x2+(34.)*x1*x1*x1*x1+(35.)*x1*x1*x1*x2+(36.)*x1*x1*x2*x2+(37.)*x1*x2*x2*x2+(38.)*x2*x2*x2*x2


0.5   1.   2.   185.375
0.5   2.   1.   596.75
0.5   3.   0.   2272.13

例子2：2021 年高教社杯全国大学生数学建模竞赛 B题 乙醇偶合制备 C4 烯烃

版权所有© Lu程序设计 2011-2021，保留所有权利
E-mail: forcal@sina.com  QQ：630715621
最近更新: 2021年11月08日