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µÈÈ¡¾ØÕóµÄÁУ¬²¢×ªÎªÒ»Î¬Êý×é
    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º¯ÊýÈ«¾ÖÓÅ»¯
};

6Â¥: Originally posted by wlfc at 2021-08-20 21:16:50
ÓÃOpenLuÇó½â£º

·ÖÎö£ºÐèÔÚ΢·Ö·½³ÌÇó½âÖд«µÝÖмä±äÁ¿lamb£¬gslÖ§³ÖµÄLuÀ©Õ¹Êýѧ¿âÖÐÇó½â΢·Ö·½³ÌµÄº¯Êýgsl_odeÌṩÁË´«µÝÖмä±äÁ¿µÄ¹¦ÄÜ¡£
½øÒ»²½ÌÖÂÛ£ºÎ¢·Ö·½³ÌÇó½â´«µÝÖмä±äÁ¿Ê±£¬°´Ê±¼ätµÄ¸Ä±äÁ¿½øÐвåÖµ ...