亚洲在线久爱草,狠狠天天香蕉网,天天搞日日干久草,伊人亚洲日本欧美

為了賬號安全,請及時綁定郵箱和手機立即綁定
已解決430363個問題,去搜搜看,總會有你想問的

Python solve_ivp 與 R lsoda

Python solve_ivp 與 R lsoda

慕無忌1623718 2021-11-16 16:28:04
我正在嘗試將 R 代碼轉換為 Python。R 代碼使用lsoda函數,它是 FORTRANDOE求解器的包裝器。Python 對應物似乎solve_ivp是 FORTRAN 的包裝器ODEPACK。我method='LSODA'在 Python 中使用它應該是 R 使用的等效方法。但是,我的結果有高達 1% 的誤差。我的代碼中沒有任何東西是隨機的,所以我相信我應該能夠完全復制結果。任何的想法?!這是 R 代碼的一部分(之前的代碼只是計算參數的值:val = c("A1" = 1, "A2" = 1, "A3" = 1, "A4" = 1, "A5" = 1, "A6" = 1, "A7" = 1)hamberg_ode <- function(t,val,p) {     dA1 = p["ktr1"]*(1 - ((p["E_MAX"] * p["C_s_gamma"])/(p["EC_50_gamma"] + p["C_s_gamma"]))) - p["ktr1"]*val["A1"]    dA2 = p["ktr1"]*val["A1"] - p["ktr1"]*val["A2"]    dA3 = p["ktr1"]*val["A2"] - p["ktr1"]*val["A3"]    dA4 = p["ktr1"]*val["A3"] - p["ktr1"]*val["A4"]    dA5 = p["ktr1"]*val["A4"] - p["ktr1"]*val["A5"]    dA6 = p["ktr1"]*val["A5"] - p["ktr1"]*val["A6"]    dA7 = p["ktr2"]*(1 - ((p["E_MAX"] * p["C_s_gamma"])/(p["EC_50_gamma"] + p["C_s_gamma"]))) - p["ktr2"]*val["A7"]    cat(val["A1"], dA1, '\n')    list(c(dA1, dA2, dA3, dA4, dA5, dA6, dA7))}out = lsoda(val, times, hamberg_ode, p)Python代碼:val = [1]*7class hamberg_ode:    def __init__(self, p):        self.p = p        return (dA1, dA2, dA3, dA4, dA5, dA6, dA7)h_function = hamberg_ode(p).fout = solve_ivp(h_function, (0, maxTime), val, t_eval=times, method='LSODA')作為數字如何發散的示例,以下是兩個代碼的 A1 和 dA1 的幾個第一個值: R1 -0.22891511 -0.22891510.9997726 -0.22879750.9997727 -0.22879760.9995454 -0.228680.9995455 -0.22868010.9901534 -0.22382210.9901523 -0.22382150.9809609 -0.21906730.9809587 -0.21906620.9719626 -0.2144130.9719604 -0.21441190.9493722 -0.20272840.9493668 -0.20272550.927996 -0.19167170.9280039 -0.19167580.9078033 -0.18122720.9078049 -0.1812280.8887056 -0.17134910.8887071 -0.1713499
查看完整描述

1 回答

?
慕俠2389804

TA貢獻1719條經驗 獲得超6個贊

正如@astoeriko 所指出的,默認相對容差 ( rtol) 在 scipy 中為 1e-3solve_ivp而在 R 中為 1e-6 lsoda。


查看完整回答
反對 回復 2021-11-16
  • 1 回答
  • 0 關注
  • 370 瀏覽
慕課專欄
更多

添加回答

舉報

0/150
提交
取消
微信客服

購課補貼
聯系客服咨詢優惠詳情

幫助反饋 APP下載

慕課網APP
您的移動學習伙伴

公眾號

掃描二維碼
關注慕課網微信公眾號