| 查看: 1446 | 回复: 5 | ||
[求助]
求助常微分方程组的参数优化问题，有没有师兄可以用1stopt帮忙试试？ 已有1人参与
|
|
师兄师姐们好！ 我最近遇到一个常微分方程组的参数优化问题。方程如下： dT <- -beta*T*V-h*F*T+rho*R dI <- beta*T*V -delta*I-k*I*F dR <- h*F*T-rho*R dV <- p*I-c*V dF <- q*I-d*F 其中： delta = 2*exp(de*(t-7)) beta,p,c,k,h,rho,q,d,de是待优化参数 初值： T = 2e7 I = 0 R = 0 V = 3 F = 1 需要拟合F和V的实验值： t V F 1 3.5 1.63 2 4.5 2.67 3 5.5 5.35 4 6.5 1.09 5 4.5 1.32 6 1.5 0 7 1.5 0 8 0 0 我尝试用r的minpack.lm包求解，但是试了两万多组初始值，始终达不到好的效果。 目前拟合最好的一组参数如下： beta =1.7e-07 p = 40 c = 17 k = 0.85 h = 0.2 rho = 1e-07 q = 3e-06 d = 0.74501237318802 de = 1 图在附件： 本人数学功底一般，刚接触这方面话题，有什么说的不清楚的地方请大家指出，我再补充。 平时潜水较多，只有10个金币啦，希望各位师兄师姐帮忙！谢谢！ 下面是我的r代码，我有说不清的大家可以看代码（解ode并画图的，不是优化的）。 #Function my_ode <- function(V_input,F_input,beta = beta, de = de, p = p, c = c, k = k, h = h, rho = rho, q = q, d = d, miu = miu, state) { library(deSolve) library(minpack.lm) library(rlist) parameters <- c(beta = beta, de = de, p = p, c = c, k = k, h = h, rho = rho, q = q, d = d, miu=miu) Target<-function(t, state, parameters) { with(as.list(c(state, parameters)),{ delta = 2*exp(de*(t-miu)) dT <- -beta*T*V-h*F*T+rho*R dI <- beta*T*V -delta*I-k*I*F dR <- h*F*T-rho*R dV <- p*I-c*V dF <- q*I-d*F list(c(dT, dI, dR, dV, dF)) }) } times <- seq(0,8,by = 0.01) out=ode(y=state,times=times,func=Target,parms=parameters) residue = c(log10(out[100,5])-3.5,log10(out[200,5])-4.5,log10(out[300,5])-5.5,log10(out[400,5])-6.5,log10(out[500,5])-4.5,log10(out[600,5])-1.5,log10(out[700,5])-1.5,out[100, 6]-1.64,out[200,6]-2.676,out[300,6]-5.35,out[400,6]-1.09,out[500,6]-1.33) residue = sqrt(sum(residue[1:7]^2)/7 + sum(residue[8:12]^2)/5) V_curve = log10(out[,5]) F_curve = out[,6] SOUTH<-1; WEST<-2; NORTH<-3; EAST<-4; GenericFigure <- function(ID, size1, size2, point, line) { plot(0:8, 0:8, type="n", xlab="X", ylab="Y"  lines(times,line, col = 2, cex = 2) points(point,col=3, pch = 15) #text(5,5, ID, col="red", cex=size1) box("plot", col="red" #mtext(paste("cex",size2,sep="" , SOUTH, line=3, adj=1.0, cex=size2, col="blue"title(paste("Figure",ID,"(",residue," ",sep=" ")} MultipleFigures <- function() { GenericFigure("V", 3, 1, V_input, V_curve) box("figure", lty="dotted", col="blue" GenericFigure("F", 3, 1, F_input, F_curve) box("figure", lty="dotted", col="blue" } par(mfrow=c(1,2)) MultipleFigures() } #Experimental data #Initial values state <- c(T = 2e7, I = 0, R = 0, V = 3, F = 1) #Parameters leiode<- function(beta,p,c,k,h,rho,q,d,de,miu) { my_ode(V_input,F_input, beta = beta, de = de, p = p, c = c, k = k, h = h, rho = rho, q = q, d = d, miu = miu, state) } #1 V_input = c(3.5,4.5,5.5,6.5,4.5,1.5,1.5,0,0,0) F_input = c(1.63,2.67,5.35,1.09,1.32) leiode(1.7e-07, 40, 17, 0.85, 0.2, 1e-07, 3e-06, 0.74501237318802,1,7)  rplot.png |
回复此楼» 猜你喜欢
光学工程学硕调剂信息
已经有32人回复
-
欢迎加入课题组
已经有0人回复
-
物理学I论文润色/翻译怎么收费?
已经有95人回复
-
散金币，求好运，祝面上顺利！
已经有78人回复
-
华东师范大学芯片设计徐珑真实水平咋样
已经有1人回复
-
调剂
已经有0人回复
-
【新加坡】纳米电子器件项目组有“联合培养博士生”名额
已经有0人回复
-
中德博士后2026项目5月30日截止- 国内博士后入站、德国两年全职
已经有0人回复
-
面向全球申请人的中德博士后2026项目5月30日截止- 国内博后入站德国两年全职
已经有0人回复
-
【新加坡】纳米电子器件项目组有“联合培养博士生”名额
已经有2人回复
-
重庆交大材料学院26年硕士招生复试实施细则已出。欢迎加入光子学微结构与器件课题组。
已经有0人回复
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感谢参与，应助指数 +1
快乐的小尼玛: 金币+10, ★★★★★最佳答案, 非常感谢！ 2016-02-26 22:30:05
感谢参与，应助指数 +1
快乐的小尼玛: 金币+10, ★★★★★最佳答案, 非常感谢！ 2016-02-26 22:30:05
|
1stOpt求这种微分方程拟合问题代码要简单的多，效果也好。 均方差(RMSE):0.85269322610406 残差平方和(SSE):11.6333718055 相关系数(R): 0.877968714324806 相关系数之平方(R^2): 0.770829063333154 修正R平方(Adj. R^2): 0.60992322741393 确定系数(DC): 0.771292669112812 F统计(F-Statistic): -1.4490058555454 参数 最佳估算 -------------------- ------------- beta1 -6.71463644100751E-8 h -7.08544427462871 rho 23.3711437688932 de -0.897079268067582 k -5.49659440033554 p 0.354150320382865 c -0.265869901817656 q 0.225238201093483 d -0.399871335562109  c1.jpg  c2.jpg |
2楼2016-02-26 16:12:27
3楼2016-02-26 22:31:45
4楼2016-02-26 23:38:29
dingd
铁杆木虫 (职业作家)
- 计算强帖: 4
- 应助: 1641 (讲师)
- 金币: 15037.3
- 散金: 101
- 红花: 234
- 帖子: 3410
- 在线: 1223.7小时
- 虫号: 291104
- 注册: 2006-10-28
5楼2016-03-02 17:15:43
6楼2016-03-02 22:45:13
