Расщепление вектора состояния при решении ряда од с использованием python odeint scipy

Question

Я конвертирую свой код из Matlab в Python и застрял на том, как разделить вектор состояния так, чтобы результат возвращал два решения.У меня есть вектор и одно значение для двух начальных условий, и в качестве конечного результата я ожидаю матрицу и вектор.

Я попытался объединить начальные условия (y0 = [c_pt_0, x_0]) в одном и том жев качестве решения (soln = [dfdt, dcdt]) (что показано ниже в коде).Я также попробовал аналогичный подход, который используется в 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()

Traceback: ValueError: установка элемента массива с последовательностью.

Любая помощь очень ценится.

РЕДАКТИРОВАНИЕ: ДОБАВЛЕННЫЙ КОД MATLAB

Вот код MATLAB, который работает, и я хочу преобразовать его в python, где вектор состояния разделен.У меня есть три файла (один основной, функция f и параметры).Прошу прощения за любые ошибки кодирования face palm , но я действительно ценю любые предложения даже для этого.1023 * и файл параметров: 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

Answer 1

После просмотра обоих кодов проблема заключается в том, что решатель odeint принимает только входы 1D-массива, а ваш 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:

* +1012 *