I can discretise the spatial domain of a system of PDEs using the method of lines, converting the system of PDEs to a system of ODEs (with a time derivative only).
These equations (for context they are mass and energy balances) are here: https://imgur.com/gallery/MYdl5KI
It can be seen that each equation on the right hand side has a number of position dependent variables (where j denotes position, i denotes species which is irrelevant here). I have seen many examples of people solving time derivative ODEs in solve_ivp/ODEint. However, the difference here is that I need to solve these equations at a series of position values also.
From what I have learned to solve for multiple positions at a single time interval I will need a for loop to set up a vector for the ODEs at each position. This vector will then be passed to solve_ivp/odeint.
My question is - how can this for loop be set up in such a way that the ODEs can be solved at a number of positions?
Any help appreciated on this
Edit: My thoughts are the vector that will be accepted by solve_ivp/odeint will look something like this:
def my_system(t, x, a, b, c, d, etc): dxdt = np.zeros(4 * n) #setting up zeros array i=0 ## Here the initial and boundary conditions are set (flow rate, temperature at inlet, initial temperature of the system, etc would be set) Tb0 = .. M0 = ... Tg0 = ... mj0 = ... y0 = ... ##The initial conditions are then used here to set up each ODE at the boundary position ##initial boundary values for the 4 ODEs dxdt[i] = ... dxdt[i + 1] = ... dxdt[i + 2] = ... dxdt[i + 3] = ... Then using the last mj equation to get the flowrate out of the boundary element mj = ... ##Solution array y will look like [y0,M0,Tg,Tb,...yn,Mn,Tgn,Tbn] at all values of t for i in range(2, 4 * n, 4): ##starting at index 2, and incrementing by 4 because we are solving for 4 values at each point dxdt[i] = ... dxdt[i + 1] = ... dxdt[i + 2] = ... dxdt[i + 3] = ... return dxdt