【答案】应助回帖
用OpenLu求解:
分析:需在微分方程求解中传递中间变量lamb,gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode提供了传递中间变量的功能。
进一步讨论:微分方程求解传递中间变量时,按时间t的改变量进行插值,可提高求解精度。例如,进行线性插值:lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)]。
另外,微分方程初值lambv未知,将其追加为拟合变量,因有5组数据,故追加5个初值lambv1,lambv2,lambv3,lambv4,lambv5。
Lu代码:CODE: !!!using["luopt","math"]; //使用命名空间
f( t,lambv,dlambv, i
: lamb
: t_A, t_lamb, t_lambdif, fac1,fac2,fac3,fac4,fac5,fac6,fac7) =
{
lamb=t_lamb[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], //线性插值
dlambv=which(t_lambdif(i)<=0.0,(lamb/3/fac7)*(fac1*((lamb/lambv)^(fac2-1)-(lambv/lamb)^(0.5*fac2+1))+fac3*((lamb/lambv)^(fac4-1)-(lambv/lamb)^(0.5*fac4+1))+fac5*((lamb/lambv)^(fac6-1)-(lambv/lamb)^(0.5*fac6+1))),(lamb/(0.001*3*fac7))*(fac1*((lamb/lambv)^(fac2-1)-(lambv/lamb)^(0.5*fac2+1))+fac3*((lamb/lambv)^(fac4-1)-(lambv/lamb)^(0.5*fac4+1))+fac5*((lamb/lambv)^(fac6-1)-(lambv/lamb)^(0.5*fac6+1)))),
0 //必须返回0
};
g( tArray1, lambv1
: lamb,lambv,max, tArray, t_pb, i,s,tf
: t_A, t_lamb, t_lambdif, fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
{
tArray=tArray1, len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_lambdif=tArray(all:2).reshape(), t_pb=tArray(all:3).reshape(), //t_A等取矩阵的列,并转为一维数组
tf=gsl_ode[@f, nil, nil, t_A, ra1(lambv1), 1e-6, 1e-6, gsl_rkf45, 1e-6,50],
i=-1, s=0.0, while{++i<max,
lamb=t_lamb(i), lambv=tf(i,1),
s=s+[fac1*(lamb^(fac2-1)*lambv^(-fac2)-lamb^(-0.5*fac2-1)*lambv^(0.5*fac2))+fac3*(lamb^(fac4-1)*lambv^(-fac4)-lamb^(-0.5*fac4-1)*lambv^(0.5*fac4))+fac5*(lamb^(fac6-1)*lambv^(-fac6)-lamb^(-0.5*fac6-1)*lambv^(0.5*fac6)) - t_pb(i)]^2.0
},
s
};
目标函数(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,lambv1,lambv2,lambv3,lambv4,lambv5
:
: fac1,fac2,fac3,fac4,fac5,fac6,fac7, tArray1,tArray2,tArray3,tArray4,tArray5)=
{
fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, //传递优化变量
g(tArray1, lambv1) + g(tArray2, lambv2) + g(tArray3, lambv3) + g(tArray4, lambv4) + g(tArray5, lambv5)
};
main(:: tArray1,tArray2,tArray3,tArray4,tArray5)=
{
tArray1=matrix{ //存放实验数据//t,lamb,lambdif,pb //0.1Hz
"0.00 0.960000 0.012566 0.057751
0.20 0.962433 0.012467 0.063943
... ... //数据省略
10.00 0.960000 0.012566 0.057751"
},
tArray2=matrix{ //存放实验数据//t,lamb,lambdif,pb //0.5Hz
"0.00 0.960000 0.062831853 0.057592
0.04 0.962137 0.062336405 0.062721
... ... //数据省略
2.00 0.960000 0.062831853 0.057592
"
},
tArray3=matrix{ //存放实验数据//t,lamb,lambdif,pb //1Hz
"0.00 0.959110 0.125664 0.101944
0.02 0.961607 0.124673 0.101015
... ... //数据省略
1.00 0.959110 0.125664 0.101944"
},
tArray4=matrix{ //存放实验数据//t,lamb,lambdif,pb //2Hz
"0.00 0.957683 0.251327 0.103554
0.02 0.962873 0.243431 0.099815
... ... //数据省略
0.50 0.957683 0.251327 0.103554"
},
tArray5=matrix{ //存放实验数据//t,lamb,lambdif,pb //5Hz
"0.00 0.953477 0.628319 0.135643
0.02 0.966487 0.508320 0.098394
... ... //数据省略
0.20 0.953477 0.628319 0.135643"
},
Opt1[@目标函数,optmax,200,optmode,10,optwaycom, optwaysimdeep, optwayconfra, optrange : 1e-5,1e10; 1e-5,1e10; 1e-5,1e10; 1e-5,1e10; -1e10,-1e-5; -1e10,-1e-5; 1e-5,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10; -1e10,1e10 ] //Opt1函数全局优化
};
结果:
4.559193187270802e-004 0.1385329162525733 2.105891636951903 1.069500926752445 -0.1172195254220325 -2.229532073322357e-004 61.97780044873748 0.943941790941042 0.9369235375692753 0.9210261290423717 0.9321935108690476 0.9479390743267424 0.4826192722906011
讨论:模型似乎与数据不匹配,拟合效果不好,且lambdif<=0与lambdif>0时公式不一样,图形上有明显转折。另外,数据稍有变动,微分方程求解结果就相差很大,例如,如果中间变量lamb未作线性插值,直接lamb=t_lamb,拟合得最优解为:
19.08839920909583 3.165001986998071e-005 0.4606063592598733 5.102787301134413 -2.111635805852794e-005 -3.981718621914052 79.73597242036105 0.9448641243550566 0.9381929974925054 0.9260059193870247 0.9399688967878371 0.9552532765459085 0.4726755802116919
两个解与变量lamb的用法互换使用时,微分方程竟然无法求解,似乎本例微分方程对数据相当敏感。
3楼求解的数据,因不知道1stopt是如何在微分方程求解中使用变量lamb的,故不知如何验证。如果将数据带入以上Lu代码求目标函数值,是非常大的。估计反过来也一样。
故微分方程求解方法不一样,获得的最优值是不一样的。对本例来说,即便变量lamb进行线性插值或未作插值,也有着天壤之别。