圖形參考：https://blog.csdn.net/wlfc/article/details/119539076

使用 3樓、4樓 的公式進行求解。

有兩個難點：
1、微分方程求解時要傳遞中間變量lamb, pa。
2、lambv為未知中間變量，需要求解。
第一個難點對Lu腳本來說不是問題，gsl支持的Lu擴展數學庫中求解微分方程的函數gsl_ode中提供了該功能。
進一步討論：微分方程求解傳遞中間變量時，按時間t的改變量進行插值，可提高求解精度。例如，進行線性插值：lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)], pa=pa0+(pa1-pa0)*[(t-t0)/(t1-t0)]。
第二個難點若使用6樓的辦法，當數據量大時似乎難以實現。故在這里使用冪級數逼近lambv，進行擬合的方法，需增加k0,k1,k2,k3四個擬合參數；當然擬合參數多時精度高，但耗時長。

仍使用部分數據擬合，以減少時間。

Lu腳本代碼：

!!!using["luopt","math"]; //使用命名空間
f(t,p,dp, i : lamb,pa, A,w,lambdif,lambvdif,lambv: t_A, t_lamb, t_pa, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4) =
{
    lamb=t_lamb+(t_lamb[i+1]-t_lamb)*[(t-t_A)/(t_A[i+1]-t_A)], pa=t_pa+(t_pa[i+1]-t_pa)*[(t-t_A)/(t_A[i+1]-t_A)], //線性插值計算lamb,pa
    lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4, //lambv=f(lamb,k0,k1,k2,k3,k4)，使用冪級數逼近
    A=0.02,w=2*pi*0.1,
    lambdif=A*w*cos(w*t),
    lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),

    dp=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
        +(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
        +fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
        -0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
        -fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
        -0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif),

    0 //必須返回0
};
目標函數(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,_k0,_k1,_k2,_k3,_k4 : i,s,tf: t_A, t_p, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
    fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, k0=_k0, k1=_k1, k2=_k2, k3=_k3, k4=_k4, //傳遞優(yōu)化變量
    //最后一個參數50表示gsl_ode函數在計算時,最多循環(huán)計算50次,這樣可以提高速度
    tf=gsl_ode[@f, nil, nil, t_A, ra1(-0.060966222), 1e-6, 1e-6, gsl_rkf45, 1e-6,50],
    sum{[tf(all:1).reshape()-t_p].^2.0} //代碼矢量化加速，特別適合大數據量
};
main(: tArray : t_A, t_lamb, t_pa, t_p)=
{
    tArray=matrix{ //存放實驗數據//t,lamb,pa,p
    "0 0.959926667 -0.218715994 -0.060966222
    0.25 0.963063333 -0.200840559 -0.046836278
    0.5 0.96612 -0.183553472 -0.035042256
    0.75 0.969016667 -0.167289626 -0.024438988
    1 0.97171 -0.152268899 -0.01706765
    1.25 0.97409 -0.139075651 -0.010939712
    1.5 0.976123333 -0.127862697 -0.006795087
    1.75 0.977786667 -0.118729759 -0.002919842
    2 0.979000001 -0.112089924 -0.001981819
    2.25 0.979743333 -0.108031322 -0.001163676
    2.5 0.98001 -0.106577018 -0.003491456
    2.75 0.979766667 -0.107904033 -0.009723697"
    },
    t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(), //t_A等取矩陣的列，并轉為一維數組
    Opt1[@目標函數, optwaysimdeep, optwayconfra] //Opt1函數全局優(yōu)化
};

結果（fac1～fac7，k0,k1,k2,k3,k4，最小值）：

17.38947868627771 0.3416680168065553 -31901.02726305072 1.871542771429994e-004 2.297270203567084e-003 8.207606707403809 -41.52220237484659 -1.705988163594953 -0.3001058349255158 2.1410790976942 4.419612473682794 -4.029380465220802 2.405714024195563e-006

繪圖代碼：

!!!using["luopt","math","win"]; //使用命名空間
f(t,p,dp, i : lamb,pa, A,w,lambdif,lambvdif,lambv: t_A, t_lamb, t_pa, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4) =
{
    lamb=t_lamb+(t_lamb[i+1]-t_lamb)*[(t-t_A)/(t_A[i+1]-t_A)], pa=t_pa+(t_pa[i+1]-t_pa)*[(t-t_A)/(t_A[i+1]-t_A)], //線性插值計算lamb,pa
    lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4, //lambv=f(lamb,k0,k1,k2,k3,k4)，使用冪級數逼近
    A=0.02,w=2*pi*0.1,
    lambdif=A*w*cos(w*t),
    lambvdif=((1.0/3/fac7)*(fac1*((lamb/lambv)^fac2-(lambv/lamb)^(0.5*fac2))+fac3*((lamb/lambv)^fac4-(lambv/lamb)^(0.5*fac4))+fac5*((lamb/lambv)^fac6-(lambv/lamb)^(0.5*fac6)))),

    dp=pa+fac1*((fac2-1)*lamb^(fac2-2)*lambv^(-fac2)*lambdif-fac2*lamb^(fac2-1)*lambv^(-fac2-1)*lambvdif
        +(0.5*fac2+1)*lamb^(-0.5*fac2-2)*lambv^(0.5*fac2)*lambdif-0.5*fac2*lamb^(-0.5*fac2-1)*lambv^(0.5*fac2-1)*lambvdif)
        +fac3*((fac4-1)*lamb^(fac4-2)*lambv^(-fac4)*lambdif-fac4*lamb^(fac4-1)*lambv^(-fac4-1)*lambvdif+(0.5*fac4+1)*lamb^(-0.5*fac4-2)*lambv^(0.5*fac4)*lambdif
        -0.5*fac4*lamb^(-0.5*fac4-1)*lambv^(0.5*fac4-1)*lambvdif)+fac5*((fac6-1)*lamb^(fac6-2)*lambv^(-fac6)*lambdif
        -fac6*lamb^(fac6-1)*lambv^(-fac6-1)*lambvdif+(0.5*fac6+1)*lamb^(-0.5*fac6-2)*lambv^(0.5*fac6)*lambdif
        -0.5*fac6*lamb^(-0.5*fac6-1)*lambv^(0.5*fac6-1)*lambvdif),

    0 //必須返回0
};
set(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,_k0,_k1,_k2,_k3,_k4 :: fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
    fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7, k0=_k0, k1=_k1, k2=_k2, k3=_k3, k4=_k4 //傳遞優(yōu)化變量
};
init(init: tArray,t_lambv,i,tf,k,lamb : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7,k0,k1,k2,k3,k4)=
{
    tArray=matrix{ //存放實驗數據//t,lamb,pa,p
    "0 0.959926667 -0.218715994 -0.060966222
    ... ...  數據省略
    2.75 0.979766667 -0.107904033 -0.009723697"
    },
    len[tArray,0,&max], t_A=tArray(all:0).reshape(), t_lamb=tArray(all:1).reshape(), t_pa=tArray(all:2).reshape(), t_p=tArray(all:3).reshape(),

    set[17.38947868627771, 0.3416680168065553, -31901.02726305072, 1.871542771429994e-004, 2.297270203567084e-003, 8.207606707403809, -41.52220237484659, -1.705988163594953, -0.3001058349255158, 2.1410790976942, 4.419612473682794, -4.029380465220802 ],

    t_lambv=new[real_s,max],
    i=-1, while{++i<max, lamb=t_lamb, t_lambv=k0+k1*lamb+k2*lamb*lamb+k3*lamb^3+k4*lamb^4},
    tf=gsl_ode[@f, nil, nil, t_A, ra1(-0.060966222), 1e-6, 1e-6, gsl_rkf45, 1e-6,50],

    cwAttach[typeSplit], cwResizePlots(1,2,2), //左二右二分裂
    k=cwAddCurve{t_A, t_p, max, 0}, //給0子圖添加曲線
    cwSetScatter(k,0),         //設置繪制點
    cwSetDataLineSize(5, k, 0),    //設置點的大小
    cwAddCurve{t_A, tf(all:1).reshape(), max, 0}, //給0子圖添加曲線
    cwAddCurve{t_A, t_lamb, max, 1}, //給1子圖添加曲線
    cwAddCurve{t_A, t_pa, max, 2},   //給2子圖添加曲線
    cwAddCurve{t_A, t_lambv, max, 3} //給3子圖添加曲線
};
ChartWnd[@init];

圖形：
https://blog.csdn.net/wlfc/article/details/119607664