我正在将代码从 Matlab 转换为 Python，并专注于如何分割状态向量，以便结果返回两个解决方案。我有一个向量和两个初始条件的单个值，我期望最终结果是一个矩阵和一个向量。
我尝试以与解决方案 (soln = [dfdt,dcdt]) 相同的方式加入初始条件 (y0 = [c_pt_0, x_0])(如下面的代码所示)。我还尝试了一种在 matlab 中使用的类似方法，即将初始条件连接到一个数组并解压结果，但我认为问题出在维度上。
#Basic imports
import numpy as np
import pylab
import matplotlib. pyplot as plt
import scipy
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# define parameters
pi = 3.14159265
V_m = 9.09
m_V__M_Pt = 1e6/195.084
rho = 21.45
R0 = 10**(-8.19)
k_d = 10**(-13)
k_r = 10**(-5)
S = 0.314 #distribution parameter
M = 0.944 #distribution parameter
## geometry
# Finite Volume Method with equidistant elements
r_max = 30.1e-9 #maximum value
n = 301 #number of elements of FVM
dr = r_max/n #length of elements, equidistant
ini_r = np.linspace(5e-10,r_max,n+1) #boundaries of elements
mid_r = ini_r[0:n]+dr/2 #center of elements
## initial conditions
#initial distribution
x0 = 1/(S*np.sqrt(2*pi)*mid_r*1e9)*np.exp((-(np.log(mid_r*1e9)-M)**2)/(2*S**2))
c_pt_0 = 0
y0 = [x0, c_pt_0]
MN_0 = scipy.trapz(np.power(mid_r, 3)*x0,
x=mid_r) # initial mass
M_0 = 4/3*pi*rho*MN_0
def f(y, t):
r = y[0]
c_pt = y[1]
#materials balance
drdt = V_m * k_r * c_pt * np.exp(-R0/ mid_r) - V_m * k_d * np.exp(R0/ mid_r)
dmdt = 4*pi*rho*mid_r**2*drdt
dMdt = np.trapz(r*dmdt, x=mid_r)
dcdt = m_V__M_Pt*(-dMdt)/M_0
dfdt = -(np.gradient(r*drdt, dr))
soln = [dfdt, dcdt]
return soln
#------------------------------------------------------
#define timespace
time = np.linspace(0, 30, 500)
#solve ode system
sln_1 = odeint (f , y0 , time,
rtol = 1e-3, atol = 1e-5)
pylab.plot(mid_r, sln_1[1,:], color = 'r', marker = 'o')
pylab.plot(mid_r, sln_1[-1,:], color = 'b', marker = 'o')
plt.show()
回溯: ValueError:使用序列设置数组元素。
非常感谢任何帮助。
编辑:添加 MATLAB 代码
这是可以运行的 MATLAB 代码，我想将其转换为 Python，其中状态向量被分割。我有三个文件(一个主文件、f 函数和参数)。请原谅任何手掌编码错误，但即使如此，我也非常感谢任何建议。
modified_model.m:
function modified_model
% import parameters
p = cycling_parameters;
% initial conditions
c_pt_0 = 0;
y0 = [p.x0; c_pt_0];
% call integrator
options_ODE=odeset('Stats','on', 'RelTol',1e-3,'AbsTol',1e-5);
[~, y] = ode15s(@(t,y) f(t, y, p), p.time, y0, options_ODE);
%% Post processing
% split state vector
r = y(:,1:p.n);
c_Pt = y(:,p.n+1);
%% Plot results
figure
hold on;
plot(p.r_m, r(1,:));
plot(p.r_m, r(end,:));
xlabel({'size'},'FontSize',15)
ylabel({'counts'},'FontSize',15)
调频
function soln = f(~, y, p)
%split state vector
r = y(1:p.n);
c_pt = y(p.n+1);
% materials balance
drdt = p.Vm_Pt.*p.k_rdp.*c_pt.*exp(-p.R0./p.r_m) - p.Vm_Pt.*p.k_dis.*exp(p.R0./p.r_m);
dmdt = 4*pi*p.rho*p.r_m.^2.*drdt;
dMdt = trapz(p.r_m, r.*dmdt);
dcdt = p.I_V*p.m_V__M_Pt*(-dMdt)/p.M_0;
dfdt = - gradient(r.*drdt,p.dr);
soln = [dfdt; dcdt];
和参数文件:cycling_parameters.m
function p=cycling_parameters
p.M = 195.084;
p.rho = 21.45;
p.time = linspace(0, 30, 500);
p.m_V__M_Pt = 1e6/p.M;
p.Vm_Pt = 9.09;
p.R0_log = -8.1963;
p.k_dis_log = -13;
p.k_rdp_log = -11;
p.R0 = 10^(p.R0_log);
p.k_dis = 10^(p.k_dis_log);
p.k_rdp = 10^(p.k_rdp_log);
%%% geometry %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Finite Volume Method with equidistant elements
p.r_max = 10.1e-9; % [m] maximum radius of PRD
p.n = 301; % number of elements of FVM
p.dr = p.r_max/p.n; % [m] length of elements, equidistant
p.r = linspace(5e-10,p.r_max,p.n+1)'; % [m] boundaries of elements
p.r_m = p.r(1:p.n)+p.dr/2; % [m] center of elements
%log normal initial distribution
S = 0.314;
M = 0.944;
p.x0 = 1./(S.*sqrt(2.*pi).*p.r_m*1e9).*exp((-(log(p.r_m*1e9)-M).^2)./(2.*S.^2));
p.r_squared = p.r_m.^2; % [m^2] squares of the radius (center of elements)
p.r_cubed = p.r_m.^3; % [m^3] cubes of the radius (center of elements)
p.MN_0 = trapz(p.r_m, p.r_cubed.*p.x0); % Eq. 2.11 denominator
p.M_0 = 4/3*pi*p.rho*p.MN_0;
p.I_V = 1; %ionomer volume fraction in the catalyst layer
最佳答案
查看这两个代码后，问题是 odeint 求解器仅接受一维数组输入，而您的 y0 是 [int, array(300,)] 而 odeint 无法使用它。但是，您可以将 y0 合并到一维数组中，然后在要积分的函数中将其拆分以进行计算，然后重新组合为输出。这是一个工作代码:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
class P:
def __init__(self, S, M):
self.M = 195.084
self.rho = 21.45
self.m_V__M_Pt = (1*10**6)/self.M
self.Vm_Pt = 9.09
self.R0_log = -8.1963
self.k_dis_log = -13
self.k_rdp_log = -11
self.R0 = 10**(self.R0_log)
self.k_dis = 10**(self.k_dis_log)
self.k_rdp = 10**(self.k_rdp_log)
self.r_max = 10.1*10**(-9)
self.n = 301
self.dr = self.r_max / self.n
self.r = np.linspace(5*10**(-10), self.r_max, self.n)
self.r_m = self.r[0:self.n+1]+self.dr/2
self.x0 = self.compute_x0(S, M)
self.r_squared = np.power(self.r_m, 2)
self.r_cubed = np.power(self.r_m, 3)
self.MN_0 = np.trapz(self.r_m, np.multiply(self.r_cubed, self.x0))
self.M_0 = (4 / 3)* np.pi * self.rho * self.MN_0
self.I_V = 1
def compute_x0(self, S, M):
p1 = np.multiply(2, np.power(S, 2))
p2 = np.multiply(S, np.sqrt(np.multiply(2, np.pi)))
p3 = np.log(self.r_m*1*10**(9)) - M
p4 = np.multiply(p2, self.r_m*10**(9))
p5 = np.power(-p3, 2)
p6 = np.multiply(p4, np.exp(np.divide(p5,p1)))
p7 = np.divide(1, p6)
return p7
def cycling_parameters():
S = 0.314
M = 0.944
p = P(S, M)
return p
def f(y, t):
p = cycling_parameters()
c_pt = y[0]
r = np.delete(y, 0)
p1 = np.multiply(p.Vm_Pt, p.k_rdp)
p2 = np.multiply(p1, c_pt)
p3 = np.multiply(p.Vm_Pt, p.k_dis)
drdt = np.multiply(p2, np.exp(np.divide(-p.R0, p.r_m))) - np.multiply(p3, np.exp(np.divide(p.R0, p.r_m)))
dmdt = np.multiply(4*np.pi*p.rho*np.power(p.r_m, 2), drdt)
p4 = np.multiply(r, dmdt)
dMdt = np.trapz(p.r_m, p4)
dcdt = p.I_V*p.m_V__M_Pt*(-dMdt)/p.M_0
p5 = np.multiply(r, drdt)
dfdt = - np.gradient(p5,p.dr)
ans = np.insert(dfdt, 0, dcdt)
return ans
def modified_model():
p = cycling_parameters()
c_pt_0 = 0
y0 = np.insert(p.x0, 0, c_pt_0)
t = np.linspace(0, 30, 500)
ans = odeint(f, y0, t, rtol = 1e-3, atol = 1e-5)
r = ans[:, 1:p.n+1]
c_Pt = ans[:, 0]
print(r)
print(c_Pt)
plt.plot(p.r_m, r[0, :], color='r', linewidth=0.5)
plt.plot(p.r_m, r[r.shape[0]-1, :], color='b', linewidth=0.5)
plt.show()
if __name__ == '__main__':
modified_model()
Python 绘图(此脚本输出的内容):
原始 Matlab 绘图:
关于python - 使用 python odeint scipy 求解一系列 ode 时分割状态向量，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/56694976/