//若长时间没有结果，在任意可接受输入的窗口，按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最优值
