Converting a matlab ODE solver to python

Question

I am attempting to convert a MATLAB code to python but I am getting answers that are completely different. I've attempted using scipy.ode,solve_ivp, and odeint.When running the code I get values that range from 1 to 0.2 but in MATLAB they range from 30 to 70. MATLAB code:

function dydt = onegroup(t,y,tsi,rho)
%Point Reactor Kinetics equation parameters
Lambda      = 10^-4;
beta        = 0.0065;
lambda      = 0.08;
%Reactivity
rho         = interp1(tsi,rho,t);

dydt = zeros(2,1);
%One Group-Delayed Precursor
dydt(1) = -lambda*y(1)+beta*y(2);
%Power
dydt(2) = ((rho-beta)/Lambda)*y(2)+(lambda*y(1))/Lambda;
end

The input file is as follows:

%%
clear;
clc;
tsi=linspace(0,20,21);
rho=ones(1,21)*0.0025;
y0= [1; 0];
ts=[0 20];

ode23(@(t,y) onegroup(t,y,tsi,rho),ts,y0)

My python code is as follows:

from scipy.integrate import ode                            
from numpy import arange,vstack,array,sqrt,ones                   
from pylab import plot,close,show,xlabel,ylabel,title,grid,pi,clf     
from scipy.interpolate import interp1d                                                    

# Function defining derivates of the positions and velocities.                                                            
def dydt(t,y,tsi,rho):

    Lambda      = 10^-4
    beta        = 0.0065
    lambda2      = 0.08

    rho = interp1d(tsi,rho, fill_value = 'extrapolate')


    #one group delayed precursor
    dydt1 = (-lambda2*y[0] + beta*y[1])
    #power
    dydt2 = (((rho(t)-beta)/Lambda)*y[1]+(lambda2*y[0])/Lambda)
    return array([dydt1, dydt2])

'''
#Final time
x = 21
#Time steps
dt = 1
tsi=np.arange(0,x,dt)
j0 = 0
times = np.arange(0,x,dt)

dt = 1

rho=np.ones(x)*0.0025
y0= [1,0]
t0 = [0,x-1]
'''

# Initial Conditions
y0, t0 = [1.,0.], 0.

# Model parameters 
k = arange(0,21) 
m = ones(21)*0.0025      

# CREATE ODE OBJECT
i = ode(dydt) # Create an ode object and bind the rhs function.
i.set_integrator('dopri5')  # Which integrator to use. 
i.set_initial_value(y0,t0)  # The initial values
# Define parameters for the derivatives function. 
# These will be passed  to the function at each time.
i.set_f_params(k,m)         

tf = 21        # Final time
dt = 1                  # Output interval
yf=[y0]                    # List for storing the output
times = arange(t0,tf,dt)    # Times to evaluate a solution. 

# Main loop for the integration
for t in times[1:]:
    i.integrate(i.t+dt)
    yf.append(i.y)

yf = array(yf)

"My python code is as follows"—cool. So what's wrong with it? Please take the tour and read How to Ask. — Chris
– Chris, Commented Apr 20, 2020 at 15:52
From simple to complex: try a very simple equation in Matlab and Scipy. Once the result is the same try something more complex. So you can rule out where the results diverge — Joe
– Joe, Commented Apr 20, 2020 at 16:29

Lutz Lehmann · Accepted Answer · 2020-04-20 18:00:37Z

2

^ in python is a bitwise logical and.

Use

Lambda = 1e-4

for that parameter.

answered Apr 20, 2020 at 18:00

Lutz Lehmann

26.1k2 gold badges26 silver badges57 bronze badges

Sign up to request clarification or add additional context in comments.

2 Comments

TylerW14 Over a year ago

The answers aren't exactly the same but they are sufficiently close now. Thank you for easily the most simple thing I have ever overlooked for too many hours!

Lutz Lehmann Over a year ago

Yes, you are using different solvers, different orders, explicit vs. implicit. What you can still do to equalize them is to explicitly set the tolerances AbsTol and RelTol for Matlab and atol and rtol in python.

Collectives™ on Stack Overflow

Converting a matlab ODE solver to python

1 Answer 1

2 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

2 Comments

Your Answer

Sign up or log in

Post as a guest

Related