python 的新功能：移植 ODE 模型

Question

我正在尝试将我的 MATLAB 代码移植到 python，但我运行遇到了一些索引问题。我在下面粘贴的代码是我计划使用 scipy.integrate.odeint 解决的一个 ode 模型。

我运行解决的第一个问题是代码 z[i-np.arange(0,(i-1)) 的索引问题（MATLAB 中的等价物是 z(i-[1:(i-1)]) ）

对 ode 模型的 python 函数的尝试：

def myModel(z,t,k,R,R2,N):

k1=k[0]
k2=k[1]
k3=k[2]

V=1/1.5*(75*z[0]+np.dot((75+57*np.arange(1,np.size(z)-2)),z[1:-2])+114*z[-2])+18*z[-1]

dzdt=np.zeros(N+2)

dzdt[0]=( 1/V*(-2*k1*z[0]*sum(z[0:-2])              
+2*k2*z[-1]*sum(z[1:-2])) ) 

dzdt[1]=( 1/V*(-2*k1*z[1]*sum(z[0:-2])          
+k1*z[0]*z[0]                               
-k2*z[1]*z[-1]                      
-k3*z[1]*V                                      
+2*k2*z[-1]*sum(z[2:-2])              
+k2*z[-2]*z[-1]) )

for i in range(2,len(z)-3):
    dzdt[i]=( 1/V*(-2*k1*z[i]*sum(z[0:-2])          
    +k1*np.dot(z[0:i],z[-(len(z)-i+1)::-1])     
    -k2*z[i]*z[-1]*i                        
    +2*k2*z[-1]*sum(z[i+1:-2])) ) 

dzdt[-3]=( 1/V*(            
+k1*np.dot(z[0:i],z[-(len(z)-i+1)::-1])
-k2*z[i]*z[-1]*i) )

dzdt[-2]=1/V*(k3*z[1]*V-k2*z[-1]*z[-2])
dzdt[-1]=1/V*(k1*sum(z[0:-2])**2-k2*z[-1]*np.dot(np.arange(0,len(z)-2),z[0:-2])+k3*z[1]*V-k2*z[-1]*z[-2]-V*(R*z[-1]/(sum(z))-R2))


return dzdt

以及 ode 模型的原始 MATLAB 函数：

function dzdt=ode(t,z,k,R,R2)


k1=k(1);
k2=k(2);
k3=k(3);


V=1/1.5*(75*z(1)+(75+57*[1:length(z)-3])*z(2:end-2)+114*z(end-1))+18*z(end);

dzdt=zeros(size(z));

dzdt(1)=1/V*(-2*k1*z(1)*sum(z(1:end-2)) ...             
    +2*k2*z(end)*sum(z(2:end-2)));                  

dzdt(2)=1/V*(-2*k1*z(2)*sum(z(1:end-2)) ...         
    +k1*z(1)*z(1) ...       
    -k2*z(2)*z(end) ...                     
    -k3*z(2)*V ...                                      
    +2*k2*z(end)*sum(z(3:end-2)) ...         
    +k2*z(end-1)*z(end));

for i=3:length(z)-3
    dzdt(i)=1/V*(-2*k1*z(i)*sum(z(1:end-2)) ...         
        +k1*z(1:i-1)'*z(i-[1:(i-1)]) ...        
        -k2*z(i)*z(end)*(i-1) ...                       
        +2*k2*z(end)*sum(z(i+1:end-2)));                       
end

dzdt(end-2)=1/V*( ...           
    +k1*z(1:i-1)'*z(i-[1:(i-1)]) ...        
    -k2*z(i)*z(end)*(i-1) ...                       
    );                          
dzdt(end-1)=1/V*(k3*z(2)*V-k2*z(end)*z(end-1));
dzdt(end)=1/V*(k1*sum(z(1:end-2))^2-k2*z(end)*[0:length(z)-3]*z(1:end-2)+k3*z(2)*V-k2*z(end)*z(end-1)-V*(R*z(end)/(sum(z))-R2));

end

提前感谢您的任何建议

Answer 1

我给你两个提示：

在Python中，您不能像在 matlab 中那样乘以向量。如果 V 是 array/list, V'*V returns matlab 中的元素乘积。在 Python (numpy) 中，您可以使用 numpy.dot 来计算两个向量的点积。
您可以使用以下格式在 numpy 中对向量进行切片：v[start:stop:step]。参见 the documentation about slicing。例如 v[::-1] 给你反向向量。

为了说明一切：

import numpy
a=numpy.array([1,2,3])
b=numpy.array([4,5,6])
a*b              # returns array([ 4, 10, 18])
numpy.dot(a,b)   # returns 32
sum(a*b)         # returns 32

结合以上提示，您不再需要该行中的 arange。

最后一点：社区很乐意帮助您解决问题，但我们不是免费编码服务。尝试展示您尝试过的内容，并尽可能少地展示您的问题，而不是给我们很多台词。这将导致更少的反对票和更多的答案。

python 的新功能：移植 ODE 模型

New to python: porting over an ODE model

python

matlab

ode