含有未知中间变量同时需要传递其他中间变量的微分方程参数拟合

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

拟合参数：fac1～fac7

未知中间变量：lambv

需传递中间变量：lamb, pa

常量：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))));

微分方程：p'=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);

数据：t,lamb,pa,p

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

分析：

1、微分方程求解时要传递中间变量：gsl支持的Lu扩展数学库中求解微分方程的函数gsl_ode中提供了该功能。

进一步讨论：微分方程求解传递中间变量时，按时间t的改变量进行插值，可提高求解精度。例如，进行线性插值：lamb=lamb0+(lamb1-lamb0)*[(t-t0)/(t1-t0)], pa=pa0+(pa1-pa0)*[(t-t0)/(t1-t0)]。

2、lambv为未知中间变量，使用幂级数逼近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[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], pa=t_pa[i]+(t_pa[i+1]-t_pa[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], //线性插值计算lamb,palambv=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, //传递优化变量//最后一个参数50表示gsl_ode函数在计算时,最多循环计算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.0609662220.25 0.963063333 -0.200840559 -0.0468362780.5 0.96612 -0.183553472 -0.0350422560.75 0.969016667 -0.167289626 -0.0244389881 0.97171 -0.152268899 -0.017067651.25 0.97409 -0.139075651 -0.0109397121.5 0.976123333 -0.127862697 -0.0067950871.75 0.977786667 -0.118729759 -0.0029198422 0.979000001 -0.112089924 -0.0019818192.25 0.979743333 -0.108031322 -0.0011636762.5 0.98001 -0.106577018 -0.0034914562.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函数全局优化
};

结果（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[i]+(t_lamb[i+1]-t_lamb[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], pa=t_pa[i]+(t_pa[i+1]-t_pa[i])*[(t-t_A[i])/(t_A[i+1]-t_A[i])], //线性插值计算lamb,palambv=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 //传递优化变量
};
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[i], t_lambv[i]=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];

图形：

含有未知中间变量同时需要传递其他中间变量的微分方程参数拟合相关推荐

python open函数encoding_python 读不同编码的文本，传递一个可选的encoding 参数给open() 函数...
文件的读写操作默认使用系统编码,可以通过调用sys.getdefaultencoding() 来得到.在大多数机器上面都是utf-8 编码.如果你已经知道你要读写的文本是其他编码方式,那么可以通过传递 ...
python之函数传递，形式参数和实际参数
python之函数传递,形式参数和实际参数. 1.函数的参数传递. 函数的参数传递有两种形式第一种是值传递,经常见于int.str.float.bool型数据.指的是将实际参数的值复制给形式参数一份 ...
python不用中间变量交换值_不使用中间变量，交换int型的 a, b两个变量的值。
不使用中间变量,交换int型的 a, b两个变量的值. 代码如下: //by ppchen var a = 10, b = 2; a = a + b; b = a - b; a = a - b; 代码 ...
（原创）matlab符号微分含有未知函数时的导数计算
我们先来看一个例子: \(y=sin(x)\),\(y\)是x的函数,同时\(x\)是关于t的函数,即为\(x(t)\),很多时候\(x(t)\)的具体表达式是未知的,这时该如何用matlab符号求\ ...
c语言函数参数传递值传递,c语言中函数参数的三种传递方式——值传递、指针传递、引用传递...
函数参数有三种传递方式值传递.指针传递.引用传递. 1.值传递将已经初始化的变量值(或常量)传递到函数中. 例如: int func(int value) { int ret = value++; ...
java 调用.net webservice axis2_java利用axis2调用.net写的webservice，传递自定义的实体类参数...
利用axis2可以很方便的自动生成客户端代码,同时对复杂参数类型的传递也很方便,本文的服务端以.net开发,有一个自定义的实体类作为参数,客户端用java,简单介绍一下利用axis2的wsdl2jav ...
java值传递和引用传递_辨析Java方法参数中的值传递和引用传递
小方法大门道小瓜瓜作为一个Java初学者,今天跟我说她想通过一个Java方法,将外部变量通过参数传递到方法中去,进行逻辑处理,方法执行完毕之后,再对修改过的变量进行判断处理,代码如下所示. publ ...
html iframe 传递数据,IFrame传入POST参数。
我有一个API,比如http://XXX/test, 只支持POST请求,需要传类似如下的参数才能拿到结果. { "id": "12", "userI ...
C++ 传递字符串数组给函数参数
C++ 传数组给一个函数,数组类型自动转换为指针类型,因而传的实际是地址. 对于传入字符串数组同理,所以如果在函数中对传入的字符串数组进行改变,函数外的字符串数组也会同时改变举个简单的例子: voi ...

含有未知中间变量同时需要传递其他中间变量的微分方程参数拟合

含有未知中间变量同时需要传递其他中间变量的微分方程参数拟合相关推荐

最新文章

热门文章