10楼: Originally posted by wlfc at 2021-08-11 16:09:33
使用 3楼、4楼 的公式进行求解。

有两个难点：
1、微分方程求解时要传递中间变量lamb, pa。
2、lambv为未知中间变量，需要求解。
第一个难点对Lu脚本来说不是问题，gsl支持的Lu扩展数学库中求解微分方程的函数 ...

7楼: Originally posted by wlfc at 2021-08-09 08:09:41
尝试用OpenLu（可从www.forcal.net下载）求解，与大家探讨。
使用 1楼 给出的第一种公式。未知中间变量lambv用解方程法求解，但不知这种方法有没有问题？
仍类似于 4楼 使用部分数据拟合，以减少运行时间。

Lu ...

!!!using["luopt","math"]; //使用命名空间

f(  t,lambv,dlambv,i
    : lamb
    : t_A, t_lamb, 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=((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)))),
    0;  //成功返回0

目标函数(_fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,lambv0
    : i,s,tf,pa,lamb,lambv
    : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
{
    fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7,  //传递优化变量
    //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
    tf=gsl_ode[@f, nil, nil, t_A, ra1(lambv0),  1e-6, 1e-6, gsl_rkf45, 1e-6,50],
    i=-1, s=0, while{++i<max,
        pa=t_pa(i), lamb=t_lamb(i), lambv=tf(i,1),
        s=s+[(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb)) - t_p(i)]^2.0
    },
    s
};

main(: tArray : t_A, t_lamb, t_pa, t_p, max)=
{
    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"
    },
    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(), //用len函数取矩阵的行数，t_A等取矩阵的列
    Opt1[@目标函数,  optwaysimdeep, optwayconfra] //Opt1函数全局优化
};

!!!using["luopt","math","win"]; //使用命名空间

f(  t,lambv,dlambv,i
    : lamb
    : t_A, t_lamb, 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=((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)))),
    0;  //成功返回0

get(t_pp,t_lambv, _fac1,_fac2,_fac3,_fac4,_fac5,_fac6,_fac7,lambv0
    : i,tf,pa,lamb,lambv
    : t_A, t_lamb, t_pa, t_p, max, fac1,fac2,fac3,fac4,fac5,fac6,fac7)=
{
    fac1=_fac1, fac2=_fac2, fac3=_fac3, fac4=_fac4, fac5=_fac5, fac6=_fac6, fac7=_fac7,  //传递优化变量
    //最后一个参数50表示gsl_ode函数在计算时,最多循环计算50次,这样可以提高速度
    tf=gsl_ode[@f, nil, nil, t_A, ra1(lambv0),  1e-6, 1e-6, gsl_rkf45, 1e-6,50],
    t_lambv.copy[tf(all,1)],
    i=-1, while{++i<max,
        pa=t_pa(i), lamb=t_lamb(i), lambv=t_lambv(i),
        t_pp[i]=(pa+fac1*((lamb/lambv)^fac2/lamb-(lambv/lamb)^(fac2*0.5)/lamb)+fac3*((lamb/lambv)^fac4/lamb-(lambv/lamb)^(fac4*0.5)/lamb)+fac5*((lamb/lambv)^fac6/lamb-(lambv/lamb)^(fac6*0.5)/lamb))
    }
};

init( main
    :  tArray, tf, i, k, t_pp,t_lambv
    : t_A, t_lamb, t_pa, t_p, max)=
{
    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"
    },
    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(), //用len函数取矩阵的行数，t_A等取矩阵的列
    t_pp=new[real_s,max], t_lambv=new[real_s,max],

    get[t_pp, t_lambv, 7.148976005697928e-003  ,  -14.5708961980305     ,    -1.059296609366114    ,    0.2342635320112431    ,    -9.68049051324744e-005  ,  -31.79361445870142     ,   -2.026215042289879     ,   1.184269268536115  ],

    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, t_pp, 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];