快樂的小尼瑪

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
帖子: 7
在線: 5.7小時
蟲號: 1948227
注冊: 2012-08-20
專業(yè): 系統(tǒng)科學(xué)與系統(tǒng)工程

[求助] 求助常微分方程組的參數(shù)優(yōu)化問題，有沒有師兄可以用1stopt幫忙試試？已有1人參與

師兄師姐們好！

  我最近遇到一個常微分方程組的參數(shù)優(yōu)化問題。方程如下：
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是待優(yōu)化參數(shù)
初值：
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包求解，但是試了兩萬多組初始值，始終達不到好的效果。
目前擬合最好的一組參數(shù)如下：
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
圖在附件：

本人數(shù)學(xué)功底一般，剛接觸這方面話題，有什么說的不清楚的地方請大家指出，我再補充。
平時潛水較多，只有10個金幣啦，希望各位師兄師姐幫忙！謝謝！
下面是我的r代碼，我有說不清的大家可以看代碼（解ode并畫圖的，不是優(yōu)化的）。
#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)

求助常微分方程組的參數(shù)優(yōu)化問題，有沒有師兄可以用1stopt幫忙試試？

rplot.png

回復(fù)此樓

» 猜你喜歡

光學(xué)工程學(xué)碩調(diào)劑信息已經(jīng)有32人回復(fù)
歡迎加入課題組已經(jīng)有0人回復(fù)
物理學(xué)I論文潤色/翻譯怎么收費? 已經(jīng)有178人回復(fù)
散金幣，求好運，祝面上順利！已經(jīng)有79人回復(fù)
華東師范大學(xué)芯片設(shè)計徐瓏真實水平咋樣已經(jīng)有1人回復(fù)
調(diào)劑已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有0人回復(fù)
中德博士后2026項目5月30日截止- 國內(nèi)博士后入站、德國兩年全職已經(jīng)有0人回復(fù)
面向全球申請人的中德博士后2026項目5月30日截止- 國內(nèi)博后入站德國兩年全職已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有2人回復(fù)
重慶交大材料學(xué)院26年碩士招生復(fù)試實施細則已出。歡迎加入光子學(xué)微結(jié)構(gòu)與器件課題組。已經(jīng)有0人回復(fù)

1樓 2016-02-26 00:32:17

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

dingd

鐵桿木蟲 (職業(yè)作家)

計算強帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時
蟲號: 291104
注冊: 2006-10-28

【答案】應(yīng)助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感謝參與，應(yīng)助指數(shù) +1
快樂的小尼瑪: 金幣+10, ★★★★★最佳答案, 非常感謝！ 2016-02-26 22:30:05

1stOpt求這種微分方程擬合問題代碼要簡單的多，效果也好。

CODE:

ConstStr delta = 2*exp(de*(tt-7));

InitialODEValue tt=0,T = 2e7,I = 0,R = 0,V = 3,F = 1;

Variable tt,V,F;

ODEFunction T' = -beta1*T*V-h*F*T+rho*R;

            I' = beta1*T*V -delta*I-k*I*F;

            R' = h*F*T-rho*R;

            V' = p*I-c*V;

            F' = q*I-d*F;

Data;

tt        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

均方差(RMSE):0.85269322610406
殘差平方和(SSE):11.6333718055
相關(guān)系數(shù)(R): 0.877968714324806
相關(guān)系數(shù)之平方(R^2): 0.770829063333154
修正R平方(Adj. R^2): 0.60992322741393
確定系數(shù)(DC): 0.771292669112812
F統(tǒng)計(F-Statistic): -1.4490058555454

參數(shù) 最佳估算
-------------------- -------------
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

求助常微分方程組的參數(shù)優(yōu)化問題，有沒有師兄可以用1stopt幫忙試試？-1

c1.jpg

求助常微分方程組的參數(shù)優(yōu)化問題，有沒有師兄可以用1stopt幫忙試試？-2

c2.jpg

贊一下

回復(fù)此樓

2樓2016-02-26 16:12:27

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

快樂的小尼瑪

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
帖子: 7
在線: 5.7小時
蟲號: 1948227
注冊: 2012-08-20
專業(yè): 系統(tǒng)科學(xué)與系統(tǒng)工程

引用回帖:

2樓: Originally posted by dingd at 2016-02-26 16:12:27
1stOpt求這種微分方程擬合問題代碼要簡單的多，效果也好。

ConstStr delta = 2*exp(de*(tt-7));
InitialODEValue tt=0,T = 2e7,I = 0,R = 0,V = 3,F = 1;
Variable tt,V,F;
ODEFunction T' = -beta1*T*V-h*F* ...

謝謝師兄！

回復(fù)此樓

3樓2016-02-26 22:31:45

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

快樂的小尼瑪

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
帖子: 7
在線: 5.7小時
蟲號: 1948227
注冊: 2012-08-20
專業(yè): 系統(tǒng)科學(xué)與系統(tǒng)工程

引用回帖:

師兄抱歉，我發(fā)帖時候遺漏了對v的描述。
v的實際值我寫錯了，應(yīng)該是：
tt       V       F
1       1e3.5       1.63
2       1e4.5       2.67
3       1e5.5       5.35
4       1e6.5       1.09
5       1e4.5       1.32
6       1e1.5       0
7       1e1.5       0
8       0       0
其中v的初始值也是不確定的，我隨便選了個3.
我畫線時候把v取了10的對數(shù)。抱歉師兄，麻煩您能幫我再試試嗎？

贊一下

回復(fù)此樓

4樓2016-02-26 23:38:29

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

dingd

鐵桿木蟲 (職業(yè)作家)

計算強帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時
蟲號: 291104
注冊: 2006-10-28

1e3.5
這是什么數(shù)值表示法？是指10^3.5？

贊一下

回復(fù)此樓

5樓2016-03-02 17:15:43

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

快樂的小尼瑪

新蟲 (初入文壇)

應(yīng)助: 0 (幼兒園)
帖子: 7
在線: 5.7小時
蟲號: 1948227
注冊: 2012-08-20
專業(yè): 系統(tǒng)科學(xué)與系統(tǒng)工程

引用回帖:

5樓: Originally posted by dingd at 2016-03-02 17:15:43
1e3.5
這是什么數(shù)值表示法？是指10^3.5？

是的師兄，平時用r比較多，r里面是這么表示的

贊一下

回復(fù)此樓

6樓2016-03-02 22:45:13

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

最具人氣熱帖推薦 [查看全部]	作者	回/看	最后發(fā)表

[考研] 材料與化工考研調(diào)劑 +9	孅華 2026-03-22	9/450	2026-03-25 13:09 by cmz0325
[考研] 303求調(diào)劑 +6	元夕元 2026-03-20	7/350	2026-03-25 12:00 by edmund7
[考研] 294求調(diào)劑材料與化工專碩 +17	陌の森林 2026-03-18	17/850	2026-03-25 11:26 by greychen00
[考研] 289求調(diào)劑 +9	懷瑾握瑜l 2026-03-20	9/450	2026-03-25 11:02 by userper
[考研] 278求調(diào)劑 +4	我可以上岸的對?/a> 2026-03-19	4/200	2026-03-25 11:01 by userper
[考研] 286求調(diào)劑 +11	Faune 2026-03-21	11/550	2026-03-25 10:11 by 霧散后相遇lc
[考研] 0854電子信息求調(diào)劑 +7	α____ 2026-03-22	8/400	2026-03-25 08:02 by 嚕gkvc
[考研] 材料學(xué)碩333求調(diào)劑 +3	北道巷 2026-03-24	3/150	2026-03-24 19:17 by pswait
[考研] 材料考研調(diào)劑生 +3	黃粱一夢千年 2026-03-24	3/150	2026-03-24 17:00 by barlinike
[考研] 求調(diào)劑 +7	十三加油 2026-03-21	7/350	2026-03-23 23:48 by 熱情沙漠
[考研] 一志愿國科過程所081700，274求調(diào)劑 +3	三水研0水立方 2026-03-23	3/150	2026-03-23 23:11 by MajorWen
[考研] 求老師收我 +3	zzh16938784 2026-03-23	3/150	2026-03-23 12:56 by ztnimte
[考研] 276求調(diào)劑 +3	YNRYG 2026-03-21	4/200	2026-03-23 08:31 by 醉在風(fēng)里
[考研] 293求調(diào)劑 +3	濤濤Wjt 2026-03-22	5/250	2026-03-22 22:21 by jiangpengfei
[考研] 生物學(xué)調(diào)劑 +5	Surekei 2026-03-21	5/250	2026-03-22 14:39 by tcx007
[考研] 材料求調(diào)劑 +5	@taotao 2026-03-21	5/250	2026-03-21 20:55 by lbsjt
[考研] 材料與化工（0856）304求B區(qū)調(diào)劑 +3	邱gl 2026-03-20	7/350	2026-03-21 19:05 by 15709483992
[考研] 考研調(diào)劑求學(xué)校推薦 +3	伯樂29 2026-03-18	5/250	2026-03-20 22:59 by JourneyLucky
[考研] 材料學(xué)碩297已過四六級求調(diào)劑推薦 +11	adaie 2026-03-19	11/550	2026-03-20 21:30 by laoshidan
[考研] 353求調(diào)劑 +3	拉鉤不許變 2026-03-20	3/150	2026-03-20 19:56 by JourneyLucky