Преобразование Matlab ODE Solver в python - PullRequest
Новая
       10

Преобразование Matlab ODE Solver в python

1 голос
/ 20 апреля 2020

Я пытаюсь преобразовать код MATLAB в python, но я получаю ответы, которые совершенно разные. Я пытался использовать scipy.ode, solve_ivp и odeint. При запуске кода я получаю значения в диапазоне от 1 до 0,2, но в MATLAB они варьируются от 30 до 70. Код MATLAB: 

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

Входной файл выглядит следующим образом: 

%%
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)

Мой python код выглядит следующим образом: 

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)

1 Ответ

1 голос
/ 20 апреля 2020

^ в python является побитовым logical and.

Используйте 

Lambda = 1e-4

для этого параметра.

...