!using["XSLSF"];                //使用命名空间XSLSF

//函数定义，用于计算微分方程组中各方程右端函数值，连分式法对微分方程组积分一步函数pbs1将调用该函数f。

//t为自变量，a1,a2,a3,a4,a6为函数值，da1,da2,da3,da4,da6为右端函数值(即微分方程的值)。

//该函数被调用时将用到参数k1,k2,k3,k4,k5,k7,k8,k9，这8个参数正好是拟合参数，用模块变量传递这8个参数。

f(t,a1,a2,a3,a4,a6,da1,da2,da3,da4,da6:x:k1,k2,k3,k4,k5,k7,k8,k9)=

{

da1=a1-k1*a1^2+k2*a3-k5*a1*a2-k7*a4*a1+k8*a6-k9*a4*a1,

da3=k1*a1^2-k2*a3-k3*a3+k4*a2,

da2=k3*a3-k4*a2-k5*a2*a1,

da4=k5*a2*a1-k7*a4*a1+k8*a6-k9*a4*a1^2,

da6=k7*a3*a1-k8*a6

};

//用连分式法对微分方程组进行积分，获得理论值。

//t1,t2为积分的起点和终点。

//h,s为自动变量。

//模块变量：hf为函数f的句柄，要预先获得该句柄；Array为工作数组；step为积分步长；eps为积分精度。

积分(t1,t2:h,s:hf,Array,step,eps)=

{

s=ceil[(t2-t1)/step],        //计算积分步数

h=(t2-t1)/s,                 //重新计算积分步长

{   pbs1[hf,t1,Array,h,eps], //对微分方程组积分一步

t1=t1+h                  //积分步长增加

}.until[abs(t1-t2)

};

//目标函数定义，自变量_k1,_k2,_k3,_k4,_k5,_k7,_k8,_k9为需要优化的参数，需要将这些参数传递给对应的模块变量k1,k2,k3,k4,k5,k7,k8,k9。

//模块变量：Array为工作数组；数组tArray存放实验数据；max为数组tArray的长度；k1,k2,k3,k4,k5,k7,k8,k9为优化变量。

优化(_k1,_k2,_k3,_k4,_k5,_k7,_k8,_k9:i,s,t0,tt,a10,a1t,a20,a2t,a30,a3t,a40,a4t,a60,a6t:Array,tArray,max,k1,k2,k3,k4,k5,k7,k8,k9)=

{

k1=_k1,k2=_k2,k3=_k3,k4=_k4,k5=_k5,k7=_k7,k8=_k8,k9=_k9, //传递优化变量，函数f中要用到k1,k2,k3,k4,k5,k7,k8,k9

i=0,s=0,

(i

tArray.getra(i*12,&t0,&tt,&a10,&a1t,&a20,&a2t,&a30,&a3t,&a40,&a4t,&a60,&a6t),

Array.setra(0,a10,a20,a30,a40,a60),   //设置积分初值，通过模块变量Array传递，Array是一个数组

积分(&t0,tt),                    //从t1积分到t2

Array.getra(0,&a10,&a20,&a30,&a40,&a60),             //从数组Array获得t2时的积分值

s=s+(a10-a1t)^2+(a20-a2t)^2+(a30-a3t)^2+(a40-a4t)^2+(a60-a6t)^2,  //累加理论值与实验值差的平方

i++

},

s

};

验证(_k1,_k2,_k3,_k4,_k5,_k7,_k8,_k9:i,s,t0,tt,a10,a1t,a20,a2t,a30,a3t,a40,a4t,a60,a6t:Array,tArray,max,k1,k2,k3,k4,k5,k7,k8,k9)=

{

k1=_k1,k2=_k2,k3=_k3,k4=_k4,k5=_k5,k7=_k7,k8=_k8,k9=_k9, //传递优化变量，函数f中要用到k1,k2,k3,k4,k5,k7,k8,k9

printff{"\r\n   tt     a1实验值      a1模拟值      a2实验值        a2模拟值      a3实验值       a3模拟值       a4实验值       a4模拟值       a6实验值       a6模拟值\r\n"},

i=0,s=0,

(i

tArray.getra(i*12,&t0,&tt,&a10,&a1t,&a20,&a2t,&a30,&a3t,&a40,&a4t,&a60,&a6t),

Array.setra(0,a10,a20,a30,a40,a60),   //设置积分初值，通过模块变量Array传递，Array是一个数组

积分(&t0,tt),                    //从t1积分到t2

Array.getra(0,&a10,&a20,&a30,&a40,&a60),             //从数组Array获得t2时的积分值

printff{"{1,r,5.2}{2,r,12.5}{3,r,18.8}{4,r,12.5}{5,r,18.8}{6,r,12.5}{7,r,18.8}{8,r,12.5}{9,r,18.8}{10,r,12.5}{11,r,18.8}\r\n",tt,a1t,a10,a2t,a20,a3t,a30,a4t,a40,a6t,a60},

i++

}

};

main(:x,xx,g,d,u,v,k,i,_eps,k1,k2,k3,k4,k5,k7,k8,k9:Array,tArray,max,hf,step,eps)=

{

hf=HFor("f"),                                    //模块变量hf保存函数f的句柄，预先用函数HFor获得该句柄

step=0.2,eps=1e-6,                               //积分步长step要合适，积分精度eps越小越精确，用于对微分方程组积分一步函数pbs1

Array=new[rtoi(real_s),rtoi(5*15)],              //申请工作数组

max=11,                                          //实验数据组数

tArray=new{rtoi(real_s),rtoi(max),rtoi(12),rtoi(EndType), //存放实验数据ti

0 , 20 , 0.8 , 0.6272 , 0 , 0.023 , 0 , 0.0159 , 0 , 0.0619 , 0 , 0.01 ,

0 , 20 , 0.6272 , 0.5769 , 0.023 , 0.0164 , 0.0114 , 0.0114 , 0.0619 , 0.1 , 0.01 , 0.02 ,

0 , 20 , 0.5769 , 0.5471 , 0.0164 , 0.0163 , 0.0103 , 0.0103 , 0.1 , 0.1213 , 0.02 , 0.026 ,

0 , 20 , 0.5471 , 0.5193 , 0.0163 , 0.0143 , 0.0092 , 0.0092 , 0.1213 , 0.1366 , 0.026 , 0.0298 ,

0 , 20 , 0.5193 , 0.4944 , 0.0143 , 0.0129 , 0.0081 , 0.0081 , 0.1366 , 0.1489 , 0.0298 , 0.0333 ,

0 , 20 , 0.4944 , 0.4673 , 0.0129 , 0.0117 , 0.0071 , 0.0071 , 0.1489 , 0.1563 , 0.0333 , 0.0287 ,

0 , 20 , 0.8 , 0.6896 , 0 , 0.0238 , 0 , 0.0138 , 0 , 0.0715 , 0 , 0.01 ,

0 , 20 , 0.6896 , 0.5716 , 0.0238 , 0.0161 , 0.0138 , 0.0104 , 0.0715 , 0.1297 , 0.02 , 0.026 ,

0 , 20 , 0.5716 , 0.5408 , 0.0161 , 0.0141 , 0.0104 , 0.0081 , 0.1297 , 0.1427 , 0.026 , 0.0298 ,

0 , 20 , 0.5408 , 0.5042 , 0.0141 , 0.0126 , 0.0081 , 0.0074 , 0.1427 , 0.1274 , 0.0298 , 0.0333 ,

0 , 20 , 0.5042 , 0.4792 , 0.0126 , 0.0116 , 0.0074 , 0.0076 , 0.1274 , 0.1625 , 0.0333 , 0.0287

},

x=new[rtoi(real_s),rtoi(9),rtoi(EndType),0.46064821368080078, 0.84128331240613474, 18.374739543487625, 4.4408607494786949, 5.2917135872283021, -1.4591430696009615, 85.986437578686022, 5.4684730589573558,0.95],  //申请工作数组，设置初值，最后一个数存放调整系数，不能为1

xx=new[rtoi(real_s),rtoi(8),rtoi(9)],            //申请工作数组

g=new[rtoi(real_s),rtoi(9)],                     //申请工作数组

_eps=1e-100, d=0,u=1.6,v=0.4,k=1000,              //变换d、u、v进一步求解，k为允许的最大迭代次数；d=0时设置的初值才有效

i=jsim[HFor("优化"),d,u,v,x,_eps,k,xx,g],        //求n维极值的单形调优法

x.getra(0,&k1,&k2,&k3,&k4,&k5,&k7,&k8,&k9,&d),   //获得最优参数值及目标终值

printff{"\r\nk={1,i}, k1={2,r}, k2={3,r}, k3={4,r}, k4={5,r}, k5={6,r}, k7={7,r}, k8={8,r}, k9={9,r}, 目标函数={10,r}\r\n",i,k1,k2,k3,k4,k5,k7,k8,k9,d},

验证(k1,k2,k3,k4,k5,k7,k8,k9),               //验证计算结果

delete[Array],delete[tArray],delete[x],delete[xx],delete[g]

};