0

I need help with solving three first-order ODEs using the scipy.integrate.ode module with the integration method of the Runge Kutta Method. My problem is that I don't know how to solve the ODE with one of the variables as a list of 1000 elements. As shown in the code section below, I want to loop through this list of 1000 values one by one when each time the ODEs are being integrated. The values "iI, e, V, Ts, Vg, Ag, N0, Tp, Gamma, Device_length" are all fixed constants. The value of "Pin_dbm" should change when each time the ODE solver is called. The unknown variables that are needed to be solved are "y[0]", "y[1]", and "dy[2]". How do I do this?

This is what I have tried:

    # Import libraries
    import matplotlib.pyplot as plt
    import numpy as np
    import os
    from scipy.integrate import ode

    # User input
    I = int(input('What is the SOA pumping current [mA]: '))        
    Device_length = float(input('What is the device (SOA) length [m]: '))

    # Fixed value constants
    iI = I*1e-3             
    Ts = 1e-9               
    Tp = 1.5e-12            
    Vg = 9e7                
    e = 1.6e-19             
    V = 50e-4*1.8e-15      
    Gamma = 0.052           
    Ag = 3e-20              
    N0 = 5e18        
    Mode  = 1       
    
    # 1000 elements list
    Pin = np.arange(1,100,0.1, dtype=float) 
    Pin_dbm = [10*np.log(P) for P in Pin]    # the list of 1000 elements
    L = np.arange(0,50,0.5, dtype=float)

    N = []                  
    S = []                  
    G = []                  
    N_end = []               
    S_end = []              
    DL = []

    def solver(Pin_dbm):
        global N
        global S
        global G

    # Function definition
        def RK4_ODE_Solver(t, y, p):

            dy = np.zeros([3])

            # three ODEs
            dy[0] = (iI/e*V) - (y[0]/Ts) - Vg*Ag*(y[0]-N0)*y[1]               
            dy[1] = y[1]/Tp + Vg*Gamma*Ag*(y[0] - N0)*y[1] + Pin_dbm          
            dy[2] = 10*np.log(np.exp(Gamma*Ag*(y[0] - N0)*Device_length))

            return dy
    
            # Initial time conditions
            t0 = 0; t_end = 1e-8; dt = 1e-13               
            y_initial = [1e16, 0]                         
            Y = [];                                
            p = [iI, e, V, Ts, Vg, Ag, N0, Tp, Gamma, [i for i in Pin_dbm]]       

            # Using RUNGE KUTTA 4th order to integrate the ODEs
            r = ode(RK4_ODE_Solver).set_integrator('dopri5', nsteps = 1e-4)
            # r = RK45(SOA_rates, t0, y_initial, t_bound=t_end-t0, first_step=1e-4)
            # r = solve_ivp(SOA_rates, method='RK45', t_span=(t0, t_end), y0=y_initial, dense_output=True) 
            r.set_f_params(p).set_initial_value(y_initial, t0)

            while r.successful() and r.t + dt < t_end:
                r.integrate(r.t + dt)
                Y.append(r.y)

            Y = np.array(Y)
            N = Y[:, 0]
            S = Y[:, 1]
            G = Y[:, 2]

            S_end.append(S[-1:])
            N_end.append(N[-1:])

    # Main function
    if (Mode == 1):
        Solver(Pin_dbm)
CYU1
  • 41
  • 5
  • Please learn about python indentation and code structure. The initialization of the stepper is never reached. The stepper itself does not perform a time loop, that you would have to do, again, see documentation. Or use the more modern steppers like `RK45` or the supplemented time-loop `solve_ivp`. Then you need some outer loop to compute each solution for each parameter separately and collect the results in one array for further processing. – Lutz Lehmann Apr 29 '22 at 07:47
  • I have tried RK45 and solve_ivp as well. None of them worked. Please be more specific about what you mean. Thanks! – CYU1 Apr 29 '22 at 15:00
  • Your use and implementation of `RK4_ODE_Solver` shows that you have problems with the basics of python. Try some easier tasks using multiple interdependent functions. – Lutz Lehmann Apr 29 '22 at 15:27
  • This code was taken reference from Nail Boohan's rate equations python code on GitHub. He coded the original code for a different purpose. I added my conditions to it to fit my purpose for the code. I followed the exact steps of a person with a PhD degree. Are you saying Nail Boohan's code doesn't follow the basic of Python? – CYU1 Apr 29 '22 at 18:17
  • Nail Boohan's Github link: https://github.com/boohann/Laser-Rate-Equations-Python/blob/master/Solving_Rate_Equations.pdf – CYU1 Apr 29 '22 at 18:37
  • The new code now indeed has the necessary time loop. The indentation, which is part of python's syntax, is still not correct. – Lutz Lehmann Apr 29 '22 at 20:35

0 Answers0