Splitting the state vector when solving a series of odes using python odeint scipy

Question

I'm converting my code from Matlab to Python and stuck on how to split the state vector such that the result returns the two solution. I have a vector and a single value for the two initial conditions and I expect as the final result a matrix and a vector.

I tried joining the initial conditions (y0 = [c_pt_0, x_0]) in the same manner as the solution (soln = [dfdt,dcdt]) (which is shown below in the code). I also tried a similar approach that is used in matlab, which is concatenating the the initial conditions to one single array and unpacking the results but I think the problem is in the dimensions.

#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: setting an array element with a sequence.

Any help is much appreciated.

EDIT: ADDED MATLAB CODE

Here is the MATLAB code that works that I want to convert to python where the state vector is split. I have three files (one main, the f function, and the parameters). Please excuse any face palm coding errors but I do appreciate any suggestions even for this.

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)

f.m

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];

and the parameters file: 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

After looking at both codes, the issue is that the odeint solver only takes 1D array inputs and your y0 is [int, array(300,)] and odeint can't work with that. However, you can merge the y0 into a 1D array and then split it up in the function you are integrating over to do the calculation then recombine as the output. Here's a working code of that:

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 plot (what this script outputs):

Original Matlab Plot: