[已解决] 微分方程组参数拟合

[复制链接]
quuq432 发表于 2016-3-1 13:34:38
本帖最后由 quuq432 于 2016-3-1 13:39 编辑

编一个关于微分方程组拟合的问题,能运行的了,也能出结果,但是得到的三个参数最佳值是错误的,且随着上下限取的不同而不同,能帮忙看看怎么办吗?
这个是我科研立项最后也是最关键一步,但是Matlab用的不熟,希望有大神帮我解决一下问题,真心感谢。
具体的微分方程是


然后求

上面L和下面L不一样,在程序里有区分。最后是得到M矩阵里M、第一个Ltsd的最佳参数值,拟合数据是t和第二个L的值(程序中已经给出了)。
代码贴出来不知道怎么回事会乱码,所以就这么直接贴了
function HH
clear all
clc
T=[502272 762325 1114354 1364122 1706100 1964980 2140916 2313059 2489034 2742752 3001044];
L=[5.65432E+42 1.30331E+43 9.00958E+42 5.3601E+42 3.56514E+42 2.5611E+42 ...
    2.27497E+42 1.94825E+42 1.50971E+42 1.22059E+42 1.04568E+42];
x0=[0.1;0.1];
M0=[0,0,0];
%M0=[0.001*1.989*10^33,1,10^42];
lb=[0.1,0.1,0.1];
ub=[0.5*1.989*10^33,1728000,1*10^46];
Lexp=...
    [502272        5.65432E+42
762325        1.30331E+43
1114354        9.00958E+42
1364122        5.3601E+42
1706100        3.56514E+42
1964980        2.5611E+42
2140916        2.27497E+42
2313059        1.94825E+42
2489034        1.50971E+42
2742752        1.22059E+42
3001044        1.04568E+42];
%t=0.1:1:402800;
[M,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@g,M0,lb,ub,[],T,x0,Lexp);

A=M(1)
B=M(2)
C=M(3)
end

function LL=g(M,T,x0,Lexp)
c=3*10^10;
k=0.5;
[t,Eint]=ode45(@f,T,x0,[],M);
LL=Eint(:,2).*t./(3.*k.*M(:,1)).*(4.*pi.*c.*Eint(:,1))-Lexp(:,2);
end

function dEint=f(t,Eint,M)   %目标函数,M(1)=M,M(2)=tsd,M(3)=L,Eint(1)代表v,Eint(2)代表Eint
c=3*10^10;
k=0.5;
dEint(1)=Eint(2)./(M(1).*Eint(1).*t);
dEint(2)=(M(3)./((1+t./M(2)).^2)-4.*pi.*c.*Eint(2).*t.*Eint(1)./(3.*k.*M(1))-Eint(2)./t);
dEint=dEint(:);
end

麻烦哪位大神能帮忙认真看一下吗。这是我科研立项最重要一步,已经卡这里很久了,真心没什么办法了。

最佳答案


fengshenone 发表于 2016-3-1 16:51:31
算了一下你给的例子,发现几个问题:
1. 你解微分方程的初值x0=[0.1;0.1],同时计算的时间是从502272开始的,换句话说你的第一个点已经确定了,那么将其带入L的表达式求,会发现无论怎么调节M的值第一个点都对不上,即L(1)不等于5.65432E+42,换句话说你的问题定义错了,或者至少你的程序错了;
2. 你这个问题其实并不好直接通过优化求解,或者说本质上是一个动态优化的问题,你这么算很容易得到局部最优解,最好对这个优化问题改进一下然后计算;
3. 个人建议可以通过调节M的值,来做Eint(:,2).*t./(3.*k.*M(:,1)).*(4.*pi.*c.*Eint(:,1))和Lexp(:,2)关于T的图,看一下各参数对整体曲线的影响,这样有助于你确定初值的选取。
4. 最后用你的方法,然后去除目标函数的第一个点,即L(1),然后从新设初值,算了一个解
1.98899999999994e+31        30893.9935465063        1.00000000000000e+46
对比了一下曲线,发现还算是那么回事,但是应该不是你想要的最优解。
P.S.多看看文献或是书类似的问题是怎么计算的,包括优化问题是怎么构建的,或是问问身边的同学。Maltab只是计算工具,算不算的出来看你对这个问题的理解和求解算法的认识
回复此楼

5 条回复


fengshenone 发表于 2016-3-1 16:51:31
算了一下你给的例子,发现几个问题:
1. 你解微分方程的初值x0=[0.1;0.1],同时计算的时间是从502272开始的,换句话说你的第一个点已经确定了,那么将其带入L的表达式求,会发现无论怎么调节M的值第一个点都对不上,即L(1)不等于5.65432E+42,换句话说你的问题定义错了,或者至少你的程序错了;
2. 你这个问题其实并不好直接通过优化求解,或者说本质上是一个动态优化的问题,你这么算很容易得到局部最优解,最好对这个优化问题改进一下然后计算;
3. 个人建议可以通过调节M的值,来做Eint(:,2).*t./(3.*k.*M(:,1)).*(4.*pi.*c.*Eint(:,1))和Lexp(:,2)关于T的图,看一下各参数对整体曲线的影响,这样有助于你确定初值的选取。
4. 最后用你的方法,然后去除目标函数的第一个点,即L(1),然后从新设初值,算了一个解
1.98899999999994e+31        30893.9935465063        1.00000000000000e+46
对比了一下曲线,发现还算是那么回事,但是应该不是你想要的最优解。
P.S.多看看文献或是书类似的问题是怎么计算的,包括优化问题是怎么构建的,或是问问身边的同学。Maltab只是计算工具,算不算的出来看你对这个问题的理解和求解算法的认识
回复此楼

quuq432 发表于 2016-3-2 22:11:49
fengshenone 发表于 2016-3-1 16:51
算了一下你给的例子,发现几个问题:
1. 你解微分方程的初值x0=[0.1;0.1],同时计算的时间是从502272开始的 ...

其实我提起结项的先交上去就是个大概算的。导师还是建议我用最小二乘法来做,让我算个精确的。我对比了网上别人的帖子类似Matlab这种编程,跟我的差不太多,但是我真的不知道我哪里错了,或者思维错了。我觉得能不能用蒙特卡洛对三个参数进行模拟,然后选一个方差最小的。

fengshenone 发表于 2016-3-3 14:43:42
quuq432 发表于 2016-3-2 22:11
其实我提起结项的先交上去就是个大概算的。导师还是建议我用最小二乘法来做,让我算个精确的。我对比了网 ...

目前来看你的模型的问题就是我说的第一点,至于你的程序本质上没问题,只是不容易得到最优解。
这个可能原因是你的参数没怎么均一化处理,数量级差的有点大,然后就是直接对微分方程这种格式进行最优化其实不太好算。当然,蒙特卡洛是肯定可以的,毕竟这类似于枚举法,而且就三个参数,网格画好一点,费点时间,肯定是能出结果的。

quuq432 发表于 2016-3-3 16:36:04
fengshenone 发表于 2016-3-3 14:43
目前来看你的模型的问题就是我说的第一点,至于你的程序本质上没问题,只是不容易得到最优解。
这个可能 ...

谢谢你啊,那我用1stopt可不可以做,写出来运行不了,你会吗?

fengshenone 发表于 2016-3-3 16:47:16
quuq432 发表于 2016-3-3 16:36
谢谢你啊,那我用1stopt可不可以做,写出来运行不了,你会吗?

没用过,不会^_^
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表