//若长时间没有结果，在任意可接受输入的窗口，按Ctrl+Alt+q退出Forcal运行。

!using["XSLSF","sys","io"]; //使用命名空间XSLSF、sys和io

//输出一维数组，该函数看不懂不要紧。将%14.6e改为%25.16e可提高输出精度

i: OutVector(p:k,i)= k=FCDLen(p),printf{"\r\n"},i=0,(i

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

//该表达式有2*n+1个参数，第一个参数为自变量，随后n个参数为函数值，最后n个参数为右端函数值(即微分方程的值)。

//n为微分方程组中方程的个数，也是未知函数的个数。

//两个冒号后的5个参数是模块变量，在这里，所有的函数属于同一个模块；模块变量在同一模块的函数内有相同的地址。

//5个模块变量k12,k21,k23,k32,k30就是拟合参数，要预先赋值，这是Forcal中通过模块变量传递参数的方法

f(t,X1,X2,X3,dX1,dX2,dX3::k12,k21,k23,k32,k30)={

dX1=-k12*X1+k21*X2,

dX2=(-k21-k23)*X2+k12*X1+k32*X3,

dX3=(-k32-k30)*X3+k23*X2

};

//用于计算目标函数：对微分方程组从t1积分到t2；x_1,x_2,x_3为实验值

//两个冒号之间的x1,x2,x3,h,i为动态变量

//两个冒号之后的都是模块变量：hf为函数f的句柄；Array为工作数组；step为积分步数；eps控制精度。这些参数要预先赋值

t_i_2(t1,t2,x_1,x_2,x_3:x1,x2,x3,h,i:hf,Array,step,eps)=

{

h=(t2-t1)/step,              //计算积分步长

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

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

}.until[abs(t1-t2)

Array.getra(0,&x1,&x2,&x3),  //从数组Array获得t2时的积分值

(x1-x_1)^2+(x2-x_2)^2+(x3-x_3)^2 //计算并返回理论值与实验值差的平方和

};

//目标函数定义。_k12,_k21,_k23,_k32,_k30为拟合参数

//两个冒号之间的t1为动态变量

//两个冒号之后的都是模块变量：Array为工作数组

//k12,k21,k23,k32,k30为拟合参数，该函数工作前，要预先赋值

J(_k12,_k21,_k23,_k32,_k30:t1:Array,k12,k21,k23,k32,k30)={

//将拟合参数_k12,_k21,_k23,_k32,_k30传给模块变量k12,k21,k23,k32,k30，在函数f中要用到k12,k21,k23,k32,k30

k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,

t1=0, Array.setra(0,100,0,0),  //设置积分初值，通过模块变量Array传递，Array是一个数组

//t1返回接近积分终值的数，理论上与积分终值相等，但实际上不一定，不过误差极小

t_i_2(&t1, 1 : 90,  8,  2)+    //从0积分到1，并计算计算理论值与实验值差的平方和

t_i_2(&t1, 3 : 73, 19,  7)+    //从1积分到3，并计算计算理论值与实验值差的平方和

t_i_2(&t1, 5 : 60, 14, 23)     //从3积分到5，并计算计算理论值与实验值差的平方和

};

//用于约束条件定义：从t1开始积分，获得t2时的积分值，由参数x1,x2,x3返回

t_i(t1,t2,x1,x2,x3:h,i:hf,Array,step,eps)=

{

h=(t2-t1)/step,

{   pbs1[hf,t1,Array,h,eps], //连分式法对微分方程组积分一步函数pbs1，hf为函数f的句柄

t1=t1+h

}.until[abs(t1-t2)

Array.getra(0,&x1,&x2,&x3)   //从数组Array获得t2时的积分值

};

//约束条件定义：n+3*m元函数，用于计算约束条件中的下限、上限及条件值：n是优化参数数目，m是约束条件数目

//_k12,_k21,_k23,_k32,_k30为优化参数；c0,c1,c2为下限值，d0,d1,d2为上限值，w0,w1,w2为条件值

S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,w0,w1,w2:t1,x1,x2,x3:Array,k12,k21,k23,k32,k30)=

{

c0=0,         c1=0,         c2=0,    //下限赋值

d0=100,       d1=100,       d2=100,  //上限赋值

//将拟合参数_k12,_k21,_k23,_k32,_k30传给模块变量k12,k21,k23,k32,k30，在函数f中要用到k12,k21,k23,k32,k30

k12=_k12,k21=_k21,k23=_k23,k32=_k32,k30=_k30,

t1=0, Array.setra(0,100,0,0),        //设置积分初值，通过模块变量Array传递，Array是一个数组

t_i(&t1, 1 :  &x1, &x2, &x3), w0=x1+x2+x3,  //条件赋值

t_i(&t1, 3 :  &x1, &x2, &x3), w1=x1+x2+x3,  //条件赋值

t_i(&t1, 5 :  &x1, &x2, &x3), w2=x1+x2+x3   //条件赋值

};

main(:a,b,alpha,_eps,k,x,xx,dx,i,tm:hf,Array,step,eps)=  //主程序，间接调用了上面定义的函数

{

tm=clock(),    //获取系统时钟脉冲

hf=HFor("f"),  //预先获得函数f的句柄，用于积分一步函数pbs1

Array=new[rtoi(real_s),rtoi(45)],  //申请工作数组，用于积分一步函数pbs1

step=20,eps=1e-6, //step为积分步数，要合适，不可太小或太大；eps越小越精确，用于积分一步函数pbs1

a=new[rtoi(real_s),rtoi(5),rtoi(EndType),1e-50,1e-50,1e-50,1e-50,1e-50], //存放常量约束条件中的变量的下界

b=new[rtoi(real_s),rtoi(5),rtoi(EndType),  1,    1,    1,    1,    1  ], //存放常量约束条件中的变量的上界

//以下数组中，前n个分量存放初始复形的第一个顶点坐标值(要求满足所有的约束条件)，返回时存放极小值点各坐标值。

//最后一个分量返回极小值。

x=new[rtoi(real_s),rtoi(6),rtoi(EndType), 0.1, 0.1, 0.1, 0.1, 0.1 ],

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

//_eps控制精度要求；alpha为反射系数；k为允许的最大迭代次数；dx为微小增量。需选择合适的alpha,_eps,k,dx值求解。

_eps=1e-10, alpha=1.01,k=800,dx=1e-4,

i=cplx[HFor("J"),HFor("S"),a,b,alpha,_eps,x,xx,k,dx], //求约束条件下n维极值的复形调优法

printf{"\r\n实际迭代次数=%d\r\n",i},

OutVector[x],    //输出数组x，最后一个数是目标函数终值，前面的是相应的优化参数值

delete[Array],delete[a],delete[b],delete[x],delete[xx] , //销毁数组

printf{"\r\n程序运行%e秒。\r\n",[clock()-tm]/1000}       //输出运行时间

};

验证初值(_k12,_k21,_k23,_k32,_k30:t1,c0,c1,c2,d0,d1,d2,w0,w1,w2:hf,Array,step,eps)=

{

hf=HFor("f"),

Array=new[rtoi(real_s),rtoi(45)],

step=20,eps=1e-6,    //step为积分步数，要合适，不可太小或太大；eps越小越精确，用于积分一步函数pbs1

S(_k12,_k21,_k23,_k32,_k30,c0,c1,c2,d0,d1,d2,&w0,&w1,&w2),

printff{"\r\n验证总量：d0={1,r}, d1={2,r}, d2={3,r}, \r\n",w0,w1,w2},

printff{"\r\n验证目标函数：J={1,r}\r\n",J(_k12,_k21,_k23,_k32,_k30)},

t1=0, Array.setra(0,100,0,0),   //设置积分初值，通过模块变量Array传递，Array是一个数组

t_i(&t1, 1 :  &c0, &c1, &c2),

printff{"\r\n实验值：t=1时，x1={1,r,-25.16} x2={2,r,-25.16} x3={3,r,-25.16} ", 90,8,2},

printff{"\r\n拟合值：t=1时，x1={1,r,-25.16} x2={2,r,-25.16} x3={3,r,-25.16} \r\n", c0,c1,c2},

t_i(&t1, 3 :  &c0, &c1, &c2),

printff{"\r\n实验值：t=1时，x1={1,r,-25.16} x2={2,r,-25.16} x3={3,r,-25.16} ", 73, 19,  7},

printff{"\r\n拟合值：t=1时，x1={1,r,-25.16} x2={2,r,-25.16} x3={3,r,-25.16} \r\n", c0,c1,c2},

t_i(&t1, 5 :  &c0, &c1, &c2),

printff{"\r\n实验值：t=1时，d0={1,r,-25.16} d1={2,r,-25.16} d2={3,r,-25.16} ", 60, 14, 23},

printff{"\r\n拟合值：t=1时，d0={1,r,-25.16} d1={2,r,-25.16} d2={3,r,-25.16} \r\n", c0,c1,c2},

delete[Array]

};

验证初值(0.1,0.1,0.1,0.1,0.1);  //第一次选取的初值

验证初值(0.112939,7.11927e-002,0.346605,1.514e-007,9.5563e-003); //forcal最优值

验证初值(0.0771532770947726, 2.59815124915171E-17, 0.197514485785414,  1.64526914364675, 3.42494200234996E-14);  //1stOpt最优值