我試圖從
uj5u.com熱心網友回復:
您必須實作系統 ODE 函式
def sys(t,u):
x,y,z1,z2 = u
v1 = z1/(p1*z2)-1
v2 = 1-y/p6
return p0*v1, p2*(x/p3-1) p4*v1, p7*z1*v2, (p7-p5*z2/z1)*z2*v2
#integrate to jump
sol1 = solve_ivp(sys,(t0,td),u0,**options)
# implement the jump
x,y,z1,z2 = sol1.y[:,-1]
fac = (2-λ)/(2 λ);
x*=fac; y*=fac
ud = [x,y,z1,z2]
# restart integration
sol2 = solve_ivp(sys,(td,tf),ud,**options)
options是方法、公差等的字典,請參閱檔案。然后評估解決方案。
uj5u.com熱心網友回復:
謝謝Lutz Lehmann,我根據您的回答解決了這個問題。附帶說明一下,等式3.3并不是多余的,我應該澄清一下,它是每枚硬幣中純銀含量的術語。所以這是一個獨立的價值(我同意在論文中定義它的奇怪方式)。
這是我最終使用的解決方案:
def dSdt(t,S):
x,y,z1,z2,z1_z2 = S
xn = p0*((z1_z2/p1)-1)
yn = p2*((x/p3)-1) p4*((z1_z2/p1)-1)
z1n = p7*z1*(1-(y/p6))
z2n = z2*(1-(y/p6))*(p7-p5*z1_z2)
z1_z2n = p5*(1-(y/p6))
return [xn,yn,z1n,z2n,z1_z2n]
# Before td
t1 = np.linspace(t0,395,abs(t0) 395,dtype=int)
S0_1 = [x_t0,y_t0,z1_t0,z2_t0,z1_div_z2_t0]
sol1 = odeint(dSdt,y0=S0_1,t=t1,tfirst=True)
# Handle division at td
sol1[-1][0] *= λ
sol1[-1][1] *= λ
S0_2 = sol1.T[:,-1]
# After td
t2 = np.linspace(start=396,stop=500,num=500-396,dtype=int)
sol2 = odeint(dSdt,y0=S0_2,t=t2,tfirst=True)
sol= np.concatenate((sol1,sol2))
t = np.concatenate((t1,t2))
轉載請註明出處,本文鏈接:https://www.uj5u.com/gongcheng/520217.html
下一篇:如何計算BST中的預期葉子數?
