快樂的小尼瑪

新蟲 (初入文壇)

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

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

師兄師姐們好！

  我最近遇到一個(gè)常微分方程組的參數(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的實(shí)驗(yàn)值：
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包求解，但是試了兩萬多組初始值，始終達(dá)不到好的效果。
目前擬合最好的一組參數(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é)功底一般，剛接觸這方面話題，有什么說的不清楚的地方請大家指出，我再補(bǔ)充。
平時(shí)潛水較多，只有10個(gè)金幣啦，希望各位師兄師姐幫忙！謝謝！
下面是我的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論文潤色/翻譯怎么收費(fèi)? 已經(jīng)有124人回復(fù)
散金幣，求好運(yùn)，祝面上順利！已經(jīng)有79人回復(fù)
華東師范大學(xué)芯片設(shè)計(jì)徐瓏真實(shí)水平咋樣已經(jīng)有1人回復(fù)
調(diào)劑已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項(xiàng)目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有0人回復(fù)
中德博士后2026項(xiàng)目5月30日截止- 國內(nèi)博士后入站、德國兩年全職已經(jīng)有0人回復(fù)
面向全球申請人的中德博士后2026項(xiàng)目5月30日截止- 國內(nèi)博后入站德國兩年全職已經(jīng)有0人回復(fù)
【新加坡】納米電子器件項(xiàng)目組有“聯(lián)合培養(yǎng)博士生”名額已經(jīng)有2人回復(fù)
重慶交大材料學(xué)院26年碩士招生復(fù)試實(shí)施細(xì)則已出。歡迎加入光子學(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è)作家)

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時(shí)
蟲號: 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)計(jì)(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小時(shí)
蟲號: 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小時(shí)
蟲號: 1948227
注冊: 2012-08-20
專業(yè): 系統(tǒng)科學(xué)與系統(tǒng)工程

引用回帖:

師兄抱歉，我發(fā)帖時(shí)候遺漏了對v的描述。
v的實(shí)際值我寫錯(cuò)了，應(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的初始值也是不確定的，我隨便選了個(gè)3.
我畫線時(shí)候把v取了10的對數(shù)。抱歉師兄，麻煩您能幫我再試試嗎？

贊一下

回復(fù)此樓

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

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

dingd

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

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時(shí)
蟲號: 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小時(shí)
蟲號: 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？

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

贊一下

回復(fù)此樓

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

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

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

[考研] 0854電子信息求調(diào)劑 +6	α____ 2026-03-22	7/350	2026-03-24 19:46 by sunjie8888
[考研] 0854人工智能方向招收調(diào)劑 +3	章小魚567 2026-03-24	3/150	2026-03-24 19:46 by zhouxuan..
[考研] 材料專碩找調(diào)劑 +5	哈哈哈吼吼吼哈 2026-03-23	5/250	2026-03-24 19:07 by 了了了了。。
[考研] 招08考數(shù)學(xué) +7	laoshidan 2026-03-20	16/800	2026-03-24 17:52 by 烏拉兒山脈
[考研] 291求調(diào)劑 +3	HanBeiNingZC 2026-03-24	3/150	2026-03-24 16:34 by barlinike
[考研] 300求調(diào)劑，材料科學(xué)英一數(shù)二 +5	leaflight 2026-03-24	5/250	2026-03-24 16:25 by laoshidan
[材料工程] 一志愿C9材料與化工專業(yè)總分300求調(diào)劑 +4	曼111 2026-03-24	5/250	2026-03-24 15:44 by 星空星月
[考研] 材料專業(yè)求調(diào)劑 +11	hanamiko 2026-03-18	11/550	2026-03-23 23:12 by peike
[考研] 一志愿陜師大生物學(xué)071000，298分，求調(diào)劑 +3	SYA！ 2026-03-23	3/150	2026-03-23 19:09 by macy2011
[考研] 求老師收我 +3	zzh16938784 2026-03-23	3/150	2026-03-23 12:56 by ztnimte
[考博] 招收博士1-2人 +3	QGZDSYS 2026-03-18	4/200	2026-03-22 10:25 by QGZDSYS
[考研] 286求調(diào)劑 +10	Faune 2026-03-21	10/500	2026-03-21 23:34 by 314126402
[考研] 材料與化工（0856）304求B區(qū)調(diào)劑 +3	邱gl 2026-03-20	7/350	2026-03-21 19:05 by 15709483992
[考研] 一志愿華中科技大學(xué)，080502，354分求調(diào)劑 +5	守候夕陽CF 2026-03-18	5/250	2026-03-21 01:06 by JourneyLucky
[考研] 304求調(diào)劑 +6	曼殊2266 2026-03-18	6/300	2026-03-21 00:32 by JourneyLucky
[考研] 296求調(diào)劑 +6	www_q 2026-03-18	10/500	2026-03-20 23:56 by JourneyLucky
[考研] 一志愿武漢理工材料工程專碩調(diào)劑 +9	Doleres 2026-03-19	9/450	2026-03-20 22:36 by JourneyLucky
[考研] 295材料求調(diào)劑，一志愿武漢理工085601專碩 +5	Charlieyq 2026-03-19	5/250	2026-03-20 20:35 by JourneyLucky
[考研] 一志愿南京航空航天大學(xué)大學(xué) ，080500材料科學(xué)與工程學(xué)碩 +5	@taotao 2026-03-20	5/250	2026-03-20 20:16 by JourneyLucky
[考研] 材料考研調(diào)劑 +3	xwt。 2026-03-19	3/150	2026-03-19 11:22 by w沐陽w