# python – System of seven ODEs solve using solve_ivp or implement RK4

I’m trying to solve a system of coupled ordinary differential equations, formed by 7 ODEs in python, using solve_ivp or either implement a fuction for RK4.

The general physical problem is as follows: Cooling of photovoltaic modules with heat exchanger coupling to the module. In this way, the module generates electrical energy and thermal energy.

I have a polynomial function, `G`

``` This function is used as input for the ODEs, which represent an electrical and a thermal system as a function of time. To solve the electrical model, I created some scripts that solve the diode equation for the photovoltaic module in question, and the output of this script is the photovoltaic power (called in the PPV thermal model) generated as a function of the module temperature and radiation . This script works great and solves part of my problem. My difficulty lies in solving the equations of the thermal model, which receives as input parameters G The equations result in this system: System of EDOS Labels: Tvidro = Tglass = T1 Tcel = Tpv = T2 Ttedlar = T3 Tabs = Tabsorber = T4 Ttubo = Ttube = T5 Tfsai = Tfluid_out = T6 Tiso = Tinsulation = T7 Using method/function for RK4, the complete code is like this (you can go direct to part "#DEFINE MODEL EQUATIONS - ODES)" : import numpy as np import matplotlib.pyplot as plt import csv from numpy.polynomial.polynomial import polyval ############################################################ with open('directory of data called teste_dados_radiacao',"r") as i: rawdata = list(csv.reader(i, delimiter = ";")) exampledata = np.array(rawdata[1:], dtype=float) xdata = exampledata[:,0] ydata = exampledata[:,1] curve = np.array(np.polyfit(xdata, ydata, 4)) rev_curve = np.array(list(reversed(curve)), dtype=float) print(rev_curve) #G_ajustado = polyval(xdata, rev_curve) """ plt.plot(xdata, ydata, label = "dados experimentais") plt.plot(xdata, model, label = "model") plt.legend() plt.show() """ ############################################################# #CONSTANTS Tamb = 25 #°C #ambient temperatura SIGMA = 5.67e-8 #W/m2K4 E_VIDRO = 0.90 #between 0.85 e 0.83 #nasrin2017 0.04 VENTO = 2 #m/s T_GROUND = Tamb + 2 #°C T_CEU = 0.00552*Tamb**1.5 Vf = 1 #m/s Do = 10e-3 #m Di = 8e-3 #m NS = 6*10 #number of cells T_F_ENT = 20 #°C #INPUTS Tcel = 25 Tv = 25 Tiso = 30 Av = 1.638*0.982 ALPHA_VIDRO = 0.9 L_VIDRO = 3e-3 #m RHO_VIDRO = 2500 #kg/m3 M_VIDRO = Av*L_VIDRO*RHO_VIDRO #kg CP_VIDRO = 500 #j/kgK K_VIDRO = 2 #W/mK TAU_VIDRO = 0.95 Pac = 0.85 H_CELL = 0.156 #m A_CELL = NS*H_CELL**2 ALPHA_CELL = 0.9 L_CEL = 3e-3 RHO_CEL = 2330 M_CEL = A_CELL*L_CEL*RHO_CEL #kg - estimated CP_CEL = 900 #J/kgK K_CEL = 140 #W/mK BETA_T = 0.43/100 # %/°C N_ELE_REF = 0.1368 #13.68% N_ELE = N_ELE_REF*(1 - BETA_T*(Tcel - 25)) #273 + 25 - tcel kelvin A_tedlar = Av L_TEDLAR = 0.33e-3 RHO_TEDLAR = 1500 M_TEDLAR = Av*L_TEDLAR*RHO_TEDLAR CP_TEDLAR = 1090 #1090 OU 2090 K_TEDLAR = 0.35 ALPHA_TEDLAR = 0.34 #doc nasa ou zero #parameters RHO_ABS = 2700 A_ABS = Av CP_ABS =900 L_ABS = 3e-3 #mm M_ABS = A_ABS*RHO_ABS*L_ABS K_ABS = 300 A_ABS_TUBO = 10*1.60*0.01+0.154*9*0.01 A_ABS_ISO = Av-A_ABS_TUBO RHO_TUBO = 2700 CP_TUBO = 900 N_TUBOS = 10 L_TUBO = N_TUBOS*1.6 M_TUBO = RHO_TUBO*L_TUBO*(3.1415/4)*(Do**2 - Di**2) K_TUBO = 300 A_TUBO_F = 0.387 #pi*Di*(L*10 VOLTAS + R(156MM)*9) A_TUBO_ISO = 0.484 #pi*Do*(L*10 VOLTAS + R(156MM)*9) A_ISO = Av RHO_ISO = 50 L_ISO = 40e-3 M_ISO = A_ISO*RHO_ISO*L_ISO CP_ISO = 670 K_ISO = 0.0375 E_ISO = 0.75 #ESTIMATED RHO_FLUIDO = 997 M_FLUIDO = L_TUBO*(3.1415/4)*Di**2*RHO_FLUIDO CP_FLUIDO = 4186 #j/kgK MI_FLUIDO = 0.890e-3 #Pa*s ou N/m2 * s K_FLUIDO = 0.607 M_PONTO = 0.05 #kg/s ou 0.5 kg/m3 #DIMENSIONLESS Pr = CP_FLUIDO*MI_FLUIDO/K_FLUIDO #water 25°C Re = RHO_FLUIDO*Vf*Di/MI_FLUIDO if (Re<=2300): Nuf = 4.364 else: Nuf = 0.023*(Re**0.8)*(Pr*0.4)*Re #COEFFICIENTS h_rad_vidro_ceu = SIGMA*E_VIDRO*(Tv**2 - T_CEU)*(Tv + T_CEU) h_conv_vidro_amb = 2.8 + 3*VENTO h_conv_tubo_fluido = 0.5*30#Nuf h_cond_vidro_cel = 1/((L_VIDRO/K_VIDRO) + (L_CEL/K_CEL)) h_cond_cel_tedlar = 1/((L_TEDLAR/K_TEDLAR) + (L_CEL/K_CEL)) h_cond_tedlar_abs = 1/((L_TEDLAR/K_TEDLAR) + (L_ABS/K_ABS)) h_cond_abs_tubo = 1/((L_TUBO/K_TUBO) + (L_ABS/K_ABS)) h_cond_abs_iso = 1/((L_ISO/K_ISO) + (L_ABS/K_ABS)) h_cond_tubo_iso = 1/((L_ISO/K_ISO) + (L_TUBO/K_TUBO)) h_conv_iso_amb = h_conv_vidro_amb h_rad_iso_ground = SIGMA*E_ISO*(Tiso**2 - T_GROUND**2)*(Tiso + T_GROUND) #GROUPS A1 = (1/(M_VIDRO*CP_VIDRO))*(ALPHA_VIDRO*Av)#*G A2 = (1/(M_VIDRO*CP_VIDRO))*(Av*(h_rad_vidro_ceu + h_conv_vidro_amb + h_cond_vidro_cel)) A3 = (1/(M_VIDRO*CP_VIDRO))*Av*h_cond_vidro_cel A4 = (1/(M_VIDRO*CP_VIDRO))*Av*(h_conv_vidro_amb + h_rad_vidro_ceu) A5 = (1/(M_CEL*CP_CEL))*(Pac*A_CELL*TAU_VIDRO*ALPHA_CELL) #*G A6 = -1*A5*N_ELE #*G A7 = (1/(M_CEL*CP_CEL))*A_CELL*h_cond_vidro_cel A8 = (1/(M_CEL*CP_CEL))*A_CELL*(h_cond_vidro_cel + h_cond_cel_tedlar) A9 = (1/(M_CEL*CP_CEL))*A_CELL*h_cond_cel_tedlar A10 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*(1 - Pac)*TAU_VIDRO*ALPHA_TEDLAR#G A11 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*(h_cond_cel_tedlar + h_cond_tedlar_abs) A12 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*h_cond_cel_tedlar A13 = (1/(M_TEDLAR*CP_TEDLAR))*A_tedlar*h_cond_tedlar_abs A14 = (1/(M_ABS*CP_ABS))*A_ABS*h_cond_tedlar_abs A15 = (1/(M_ABS*CP_ABS))*(A_ABS*h_cond_tedlar_abs + A_ABS_TUBO*h_cond_abs_tubo + A_ABS_ISO*h_cond_abs_iso) A16 = (1/(M_ABS*CP_ABS))*A_ABS_TUBO*h_cond_abs_tubo A17 = (1/(M_ABS*CP_ABS))*A_ABS_ISO*h_cond_abs_iso A18 = (1/(M_TUBO*CP_TUBO))*A_ABS_TUBO*h_cond_abs_tubo A19 = (1/(M_TUBO*CP_TUBO))*(A_ABS_TUBO*h_cond_abs_tubo + A_TUBO_F*h_conv_tubo_fluido + A_TUBO_ISO*h_cond_tubo_iso) A20 = (1/(M_TUBO*CP_TUBO))*A_TUBO_F*h_conv_tubo_fluido*0.5 A21 = (1/(M_TUBO*CP_TUBO))*A_TUBO_ISO*h_cond_tubo_iso A22 = (1/(M_FLUIDO*CP_FLUIDO))*A_TUBO_F*h_conv_tubo_fluido A23 = (1/(M_FLUIDO*CP_FLUIDO))*(A_TUBO_F*h_conv_tubo_fluido*0.5 + M_PONTO*CP_FLUIDO) A24 = (1/(M_FLUIDO*CP_FLUIDO))*(T_F_ENT*(M_PONTO*CP_FLUIDO - h_conv_tubo_fluido*A_TUBO_F*0.5)) A25 = (1/(M_ISO*CP_ISO))*A_ABS_ISO*h_cond_abs_iso A26 = (1/(M_ISO*CP_ISO))*(A_ABS_ISO*h_cond_abs_iso + A_TUBO_ISO*h_cond_tubo_iso + A_ISO*h_conv_iso_amb + A_ISO*h_rad_iso_ground) A27 = (1/(M_ISO*CP_ISO))*A_TUBO_ISO*h_cond_tubo_iso A28 = (1/(M_ISO*CP_ISO))*A_ISO*(h_conv_iso_amb*Tamb + h_rad_iso_ground*T_GROUND) #DEFINE MODEL EQUATIONS - ODES - (GLASS, PV CELL, TEDLAR, ABSORBER, TUBE, FLUID, INSULATION) # dT1dt = A1*G_ajustado - A2*x[0] + A3*x[1] + A4 # dT2dt = A5*G_ajustado - A6*G_ajustado + A7*x[0] - A8*x[1] + A9*x[2]# dT3dt = A10*G_ajustado - A11*x[2] + A12*x[1] +A13*x[3] def SysEdo(x, k):#tv-x[0] tcel-x[1] ttedlar-x[2] tabs-x[3] ttubo-x[4] tiso-x[5] tfs-x[6] dT1dt = A1*polyval(k,rev_curve) - A2*x[0] + A3*x[1] + A4 dT2dt = A5*polyval(k,rev_curve) - A6*polyval(k,rev_curve) + A7*x[0] - A8*x[1] + A9*x[2] dT3dt = A10*polyval(k,rev_curve) - A11*x[2] + A12*x[1] +A13*x[3] dT4dt = A14*x[2] - A15*x[3] + A16*x[4] + A17*x[5] dT5dt = A18*x[3] - A19*x[4] + A20*x[6] + A20*T_F_ENT + A21*x[5] dT6dt = A22*x[4] - A23*x[6] + A24 dT7dt = A25*x[3] - A26*x[5] + A27*x[4] + A28 Tdot = np.array([dT1dt, dT2dt, dT3dt, dT4dt, dT5dt, dT6dt, dT7dt]) return Tdot #RungeKutta4 def RK4(f, x0, t0, tf, dt): t = np.arange(t0, tf, dt) #time vector nt = t.size #lenght of time vector nx = x0.size #length of state variables? x = np.zeros((nx,nt)) #initialize 2D vector x[:,0] = x0 #initial conditions #RK4 constants for k in range(nt-1): k1 = dt*f(t[k], x[:,k],k) k2 = dt*f(t[k] + dt/2, x[:,k] + k1/2, k) k3 = dt*f(t[k] + dt/2, x[:,k] + k2/2, k) k4 = dt*f(t[k] + dt, x[:,k] + k3, k) dx = (k1 + 2*k2 + 2*k2 + k4)/6 x[:,k+1] = x[:,k] + dx return x,t #Define problems f = lambda t, x, k : SysEdo(x, k) #initial state - t0 is initial time - tf is final time - dt is time step x0 = np.array([30, 30, 30, 30, 30, 30, 30]) t0 = 0 tf = 1000 dt = 1 #EDO SOLVE x, t = RK4(f, x0, t0, tf, dt) plt.figure() plt.plot(t, x[0], '-', label="Tvidro") """ plt.plot(t, x[1], '-', label="Tpv") plt.plot(t, x[2], '-', label="Ttedlar") plt.plot(t, x[3], '-', label="Tabs") plt.plot(t, x[4], '-', label="Tiso") plt.plot(t, x[5], '-', label="Ttubo") plt.plot(t, x[6], '-', label="Tfsai")""" plt.title('Gráfico') plt.legend(['Tvidro', 'Tpv', 'Ttedlar', 'Tabs', 'Tiso', 'Ttubo', 'Tfsai'], shadow=False) plt.xlabel('t (s)') plt.ylabel('Temperatura (°C)') plt.xlim(0,20) plt.ylim(0,150) plt.grid('on') plt.show() Thank you in advance, I am also open to completely start the implementation from scratch if there is a better way to do this with python or matlab. ```
``` ```
``` Categories Coding Forums Tags implement, ODEs, python, RK4, solve, solveivp, System Post navigation Efficient Heterogeneous Parallel Programming Using OpenMPoneAPI – The Cross-Architecture, Multi-Vendor Path to Accelerated Computing ```
``` ```
``` Leave a Comment CommentName Email Website Save my name, email, and website in this browser for the next time I comment. ```
``` ```
``` SearchRecent Postscurl – Poetry – RuntimeError: No python executable found in shell environment. Tried: [‘python3’, ‘python’, ‘py.exe -3’, ‘py.exe -2’] Mid-priced robot vacuum and mop with auto-empty station Python operator is not being called in dynamic subdag in Airflow These are the Best Kindle E-readers to buy in 2022 The best new features in .NET 6 Recent CommentsNo comments to show. ```
``` ```
``` About UsContact UsDMCAPrivacy PolicyTerms & Conditions © 2022 pcodmiracles • Built with GeneratePress !function(){"use strict";if("querySelector"in document&&"addEventListener"in window){var e=document.body;e.addEventListener("mousedown",function(){e.classList.add("using-mouse")}),e.addEventListener("keydown",function(){e.classList.remove("using-mouse")})}}();.wp-container-1 > .alignleft { float: left; margin-inline-start: 0; margin-inline-end: 2em; }.wp-container-1 > .alignright { float: right; margin-inline-start: 2em; margin-inline-end: 0; }.wp-container-1 > .aligncenter { margin-left: auto !important; margin-right: auto !important; } .wp-container-2 > .alignleft { float: left; margin-inline-start: 0; margin-inline-end: 2em; }.wp-container-2 > .alignright { float: right; margin-inline-start: 2em; margin-inline-end: 0; }.wp-container-2 > .aligncenter { margin-left: auto !important; margin-right: auto !important; } var wpcf7 = {"api":{"root":"https:\/\/pcodmiracles.com\/wp-json\/","namespace":"contact-form-7\/v1"},"cached":"1"}; var generatepressMenu = {"toggleOpenedSubMenus":"1","openSubMenuLabel":"Open Sub-Menu","closeSubMenuLabel":"Close Sub-Menu"}; !function(){window.advanced_ads_ready_queue=window.advanced_ads_ready_queue||[],advanced_ads_ready_queue.push=window.advanced_ads_ready;for(var d=0,a=advanced_ads_ready_queue.length;d<a;d++)advanced_ads_ready(advanced_ads_ready_queue[d])}(); ```