Search a 2D Matrix II(C++搜索二维矩阵 II)

清除浮动方法梳理

  返回  

初值未知且需要传递中间变量的微分方程参数拟合

2021/8/20 20:22:36 浏览:

问题来源:1stopt 7参数拟合求助 - 计算模拟 - 程序代码 - 小木虫论坛-学术科研互动平台 (muchong.com)

根据两个式子(一个微分方程,一个代数方程)拟合7个未知参数fac1~fac7。微分方程中lamb为需要传递的中间变量,另外,初值lambv未知。

lambv'=((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))));

p=(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

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

分析:gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode提供了传递中间变量的功能。
进一步讨论:微分方程求解传递中间变量时,按时间t的改变量进行插值,可提高求解精度。例如,进行线性插值:lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)]。

另外,微分方程初值lambv未知,将其追加为拟合变量lambv0。

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函数全局优化
};

结果:

7.148976005697928e-003    -14.5708961980305         -1.059296609366114        0.2342635320112431        -9.68049051324744e-005    -31.79361445870142        -2.026215042289879        1.184269268536115         3.55377979903292e-006

绘图代码:

!!!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];

图形:

 

联系我们

如果您对我们的服务有兴趣,请及时和我们联系!

服务热线:18288888888
座机:18288888888
传真:
邮箱:888888@qq.com
地址:郑州市文化路红专路93号