Perse95 Code

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import math


def initialize_system(*constants):
    #sets the appropriate initial/boundary conditions

    dx, number_of_species , number_of_nodes, A_s, C_pj = constants #unpack constants
    x0_species_eqn = np.zeros((2,number_of_species, number_of_nodes))
    x0_other_eqn = np.zeros(3,number_of_nodes)

    #here the initial and boundary conditions are set (flow rate, Tin, T0, etc)
    Tb0 = ...
    M0 = ...
    Tg0 = ...
    mj0 = ...
    y0 = ...

    for i in range(number_of_species):
        for j in range(number_of_nodes):
            #For each species i , we set the appropriate initial state for each node j for the dyi/dt Equation
            x0_species_eqn[0][i][j] = ...
            #Similarly for the dM_ij/dt equation
            x0_species_eqn[1][i][j] = ...
            mj = mj(a,b,c) #calculate mj

    for j in range(number_of_nodes):
        # Similarly for the dT_g/dt equation
        x0_other_eqn[0][j] = ...
        # Similarly for the dT_b/dt equation
        x0_other_eqn[1][j] = ...

        mj = mj(a,b,c) #calculate mj

    return np.concatenate([x0_species_eqn.reshape(-1), x0_other_eqn.reshape(-1)], axis=0)


def df(t,x,*constants):
    #the time derivative of x
    dx, number_of_species, number_of_nodes, A_s, C_pj = constants  # unpack constants

    #we operate on the array as before, but first we unpack the initial state
    x_species_eqn = x.reshape(2,number_of_species,number_of_nodes)
    x_other_eqn = x.reshape(3,number_of_nodes)

    dxdt_species_eqn = np.zeros_like(x_species_eqn)
    dxdt_other_eqn = np.zeros_like(x_other_eqn)

    for i in range(number_of_species):
        for j in range(number_of_nodes):
            # For each species i , we compute dyi/dt Equation for each node j
            dxdt_species_eqn[0][i][j] = ...
            # Similarly for the dM_ij/dt equation
            dxdt_species_eqn[1][i][j] = ...
            # Similarly for the dT_g/dt equation

            mj = mj(a,b,c)

    for j in range(number_of_nodes):
        dxdt_other_eqn[0][j] = ...
        # Similarly for the dT_b/dt equation
        dxdt_other_eqn[1][j] = ...

        mj = mj(a,b,c)


            #what are these equations called (the relevant ODE for each species and position?)
    return np.concatenate([dxdt_species_eqn.reshape(-1), dxdt_other_eqn.reshape(-1)], axis=0)

def mj(a, b, c):
    #this function calculates mj as a function of the relevant paramteters
    #this function will be called inside the initialize_system and df functions where necessary

n = 10
A_s = 1
C_pj = 1 #can be an array

#create your initial conditions
t0 = 0
tf = 100 #how long you want to integrate for
x0 = initialize_system(n,A_s,C_pj)

#solve the system
sol = solve_ivp(df, [t0,tf], x0, dense_output = True) #dense_output gives you an interpolant (3rd/4th order polynomial)

#decide at which times you'd like to check your system
t_eval = [t_i, t_i+1,...t_f]

x_eval = [sol.sol(t) for t in t_eval]


x_species_eqn, x_other_eqn = x[:2*number_of_species*number_of_nodes], x[2*number_of_species*number_of_nodes:]

x_species_eqn = x_species_eqn.reshape(2, number_of_species, number_of_nodes)
x_other_eqn   = x_other_eqn.reshape(3, number_of_nodes)


# And then use them like you normally would
...