4樓: Originally posted by wlfc at 2021-08-18 08:07:39
lambv'=if(lambdif<=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* ...

!!!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等取矩陣的列，并轉(zhuǎn)為一維數(shù)組
    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
};
目標(biāo)函數(shù)(_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,  //傳遞優(yōu)化變量
    g(tArray1, lambv1) +  g(tArray2, lambv2) +  g(tArray3, lambv3) +  g(tArray4, lambv4) +  g(tArray5, lambv5)
};
main(:: tArray1,tArray2,tArray3,tArray4,tArray5)=
{
    tArray1=matrix{ //存放實(shí)驗(yàn)數(shù)據(jù)//t,lamb,lambdif,pb  //0.1Hz
    "0.00         0.960000         0.012566         0.057751
0.20         0.962433         0.012467         0.063943
... ...  //數(shù)據(jù)省略
10.00         0.960000         0.012566         0.057751"
    },

    tArray2=matrix{ //存放實(shí)驗(yàn)數(shù)據(jù)//t,lamb,lambdif,pb  //0.5Hz
    "0.00         0.960000         0.062831853        0.057592
0.04         0.962137         0.062336405        0.062721
... ...  //數(shù)據(jù)省略
2.00         0.960000         0.062831853        0.057592
"
    },

    tArray3=matrix{ //存放實(shí)驗(yàn)數(shù)據(jù)//t,lamb,lambdif,pb  //1Hz
    "0.00         0.959110         0.125664         0.101944
0.02         0.961607         0.124673         0.101015
... ...  //數(shù)據(jù)省略
1.00         0.959110         0.125664         0.101944"
    },

    tArray4=matrix{ //存放實(shí)驗(yàn)數(shù)據(jù)//t,lamb,lambdif,pb  //2Hz
    "0.00         0.957683         0.251327         0.103554
0.02         0.962873         0.243431         0.099815
... ...  //數(shù)據(jù)省略
0.50         0.957683         0.251327         0.103554"
    },

    tArray5=matrix{ //存放實(shí)驗(yàn)數(shù)據(jù)//t,lamb,lambdif,pb  //5Hz
    "0.00         0.953477         0.628319         0.135643
0.02         0.966487         0.508320         0.098394
... ...  //數(shù)據(jù)省略
0.20         0.953477         0.628319         0.135643"
    },

    Opt1[@目標(biāo)函數(shù),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函數(shù)全局優(yōu)化
};

6樓: Originally posted by wlfc at 2021-08-20 21:16:50
用OpenLu求解：

分析：需在微分方程求解中傳遞中間變量lamb，gsl支持的Lu擴(kuò)展數(shù)學(xué)庫(kù)中求解微分方程的函數(shù)gsl_ode提供了傳遞中間變量的功能。
進(jìn)一步討論：微分方程求解傳遞中間變量時(shí)，按時(shí)間t的改變量進(jìn)行插值 ...