用1stopt拟合常微分方程参数
求助各位大佬,帮忙用1stopt高版本跑下面的代码,感激不尽!!!
title "canshunihe";
constant yield=0.229, k=60.328, mu=2.282,b=2.854, q=3.933, n=1.4,alpha=73.408, gamma0=5.751, beta=2.188, phai=0.141,e=1000;
//initialodevalue t=[0,8000],y1=[0,], y2=[0,], y3=[8,],y4=[0.229,],y5=[0.001,];
variable t,y3,y4;
parameter c1,c2,c3;
//odeoptions = [s=1,a=4,p=50]; //s: step size = 1, a: algorithm (1: euler, 2: rk2, 3: rk3, 4: rk4), p: population size = 5;
plot y4,y5;
odefunction y1'=((y4/(1-y5)-y2-yield)/c1)^c2*y3^(-mu);
y2'=b*(q-y2)*y1';
y3'=alpha*y3^(-gamma0)+beta*y1'*y3^(-phai);
y4'=e*(1e-3-y1');
y5'=y5^c3;
// t y4
rowdata;
0,100,200,300,400,500,600,700,800,900,1000,1500,2000,3000,4000,5500,6500,7500,8000;
8,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.0,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8;
0,2.63,2.95,2.99,3.52,3.81,4.05,4.25,4.43,4.6,4.72,5.07,5.28,5.2,5.01,4.65,4.2,3.55,1.04;
返回小木虫查看更多
参考下,微分方程拟合比较费时间,不保证是最优解:
均方差(RMSE): 0.532681607689156
残差平方和(SSR): 10.2149890261309
相关系数(R): 0.932921869702266
相关系数之平方(R^2): 0.870343214968771
修正R平方(Adj. R^2): 0.845683647429779
确定系数(DC): 0.718147888636154
F统计(F-Statistic): 28.6980923726451
参数 最佳估算
-------------------- -------------
c1 -3.77485429494E-5
c2 -0.341995932934685
c3 1.1560796736969
y1 初值 4.64067042228588
y2 初值 2.77131321062938
y5 初值 2.15540575432628E-14
下面效果更好些:
均方差(RMSE): 0.42661365444081
残差平方和(SSR): 6.55197156559236
相关系数(R): 0.980897479289642
相关系数之平方(R^2): 0.962159864876774
修正R平方(Adj. R^2): 0.954431386315503
确定系数(DC): 0.778949690316621
F统计(F-Statistic): 22.6179953241421
参数 最佳估算
-------------------- -------------
c1 1.19861624726491
c2 32.5275008605798
c3 2.24973236410188
y1 初值 -1.56792343709876
y2 初值 -0.697815532746156
y5 初值 0.00313654380605316
感谢dingd大神的回复。
我想要得到的结果是考虑断裂的应力应变曲线图,类似下图所示。走势应该是先上升后下降。我用您得到的参数,在matlab计算的时候,有些不符合实际的曲线。
可否让我看下拟合的结果图,再次感谢。
结果图.jpg
还有就是,我仅仅知道y1和y2的初值为0,y5的初值为0.001.这样的话计算的值可能就会更符合实际,但是我又不知道怎么设定他们的初值。
能在帮我把下面的代码算一下么?
title "canshunihe";
constant yield=0.229, k=60.328, mu=2.282,b=2.854, q=3.933, n=1.4,alpha=73.408, gamma0=5.751, beta=2.188, phai=0.141,e=1000;
initialodevalue t=0,y1=0, y2=0, y3=8,y4=0.229,y5=0.001;
variable t,y3,y4;
parameter c1,c2,c3;
//odeoptions = [s=1,a=4,p=50]; //s: step size = 1, a: algorithm (1: euler, 2: rk2, 3: rk3, 4: rk4), p: population size = 5;
plot y4,y5;
odefunction y1'=((y4/(1-y5)-y2-yield)/c1)^c2*y3^(-mu);
y2'=b*(q-y2)*y1';
y3'=alpha*y3^(-gamma0)+beta*y1'*y3^(-phai);
y4'=e*(1e-3-y1');
y5'=y5^c3;
// t y4
rowdata;
t=0,100,200,300,400,500,600,700,800,900,1000,1500,2000,3000,4000,5500,6500,7500,8000;
y3=8,8.1,8.2,8.3,8.4,8.5,8.6,8.7,8.8,8.9,9.0,9.1,9.2,9.3,9.4,9.5,9.6,9.7,9.8;
y4=0.229,2.63,2.95,2.99,3.52,3.81,4.05,4.25,4.43,4.6,4.72,5.07,5.28,5.2,5.01,4.65,4.2,3.55,1.04;
Parameters k1,k2,k3,k4,k5,k6,k7;
InitialODEValue t=0,a=0,b=0,c=0,d=0,e=0.0000074;
Variable t,a,b,c,d,e;
ODEFunction e'=-k1*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5-k2*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5-k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5;
a'=6*k1*3051*e*172/0.000001343*(1-1.13*Ln(e/0.0000074))^0.5+2*k2*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k4*b*3051-k5*a*d-3*k6*a*b+3*k7*c*3051;
b'=4*k2*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+2*k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343-k4*b*3051+k5*a*d-k6*a*b+k7*c*3051;
c'=k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k6*a*b-k7*c*3051;
d'=4*k1*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k3*3051*e*172*(1-1.13*Ln(e/0.0000074))^0.5/0.000001343+k4*3051*b-k5*d*a;
Data;
//t H2 a CO b CH4 c CO2 d
60 5.317534543 2.282574774 0.076404295 9.720719792
300 16.53258605 0.921333731 0.460228388 14.65292426
600 21.35578457 0.748749046 0.707176158 17.85455396
1800 30.09578114 0.572388504 1.333230207 22.92109313
2400 32.48226677 0.445459631 1.986111389 26.69658144
有没有好心人帮忙运行一下这个程序呀,谢谢啦~,