How can I fix ODE45 returning a straight line when using specific `tspan` points, instead of showing dynamic behavior over the time range?

Question

Ehtisham 2024-9-20，9:16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2154185-how-can-i-fix-ode45-returning-a-straight-line-when-using-specific-tspan-points-instead-of-showing

评论： Sam Chak 2024-9-23，11:33

6 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2024-9-20，15:46

编辑：Torsten 2024-9-20，15:49

在 MATLAB Online 中打开

The code below works, but I doubt it is really what you want. Try to explain yourself more clearly (maybe by using the Google Translator).

% Define parameters

Kf_Max = 100; % maximum forward rate

Kb_Max = 80; % maximum backward rate

Kb_Min = 10 ;

RLC = [0.1, 0.5, 10, 1]; % RLC values

TauKf_ON = -0.1; % TauKf_ON

TauKf_OFF = -0.1; % TauKf_OFF

TauKb_ON = -0.01; % TauKb_ON

TauKb_OFF = -0.01; % TauKb_ON

PhaseTimes = [0,10, 1000, 1500, 2000];% phase times (each row defines a phase)

% Example initial conditions for non-active and active receptor concentrations

Non_Active_Receptor_concentration = 100; % Initial non-active receptor concentration (as a vector)

Active_Receptor_concentration = 0; % Initial active receptor concentration (as a vector)

[t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase (Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...

Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes,Non_Active_Receptor_concentration, Active_Receptor_concentration);

function [t, Non_Active_Receptor_concentration, Active_Receptor_concentration, PhaseResults, Kf_LMax, Kb] = SimfileNPhase(Kf_Max, RLC, TauKf_ON, TauKf_OFF, ...

Kb_Max, Kb_Min, TauKb_ON, TauKb_OFF, PhaseTimes, Non_Active_Receptor_concentration, Active_Receptor_concentration)

% Ensure that Non_Active_Receptor_concentration and Active_Receptor_concentration are column vectors

Non_Active_Receptor_concentration = Non_Active_Receptor_concentration(:);

Active_Receptor_concentration = Active_Receptor_concentration(:);

% Set initial conditions

initial_conditions = [Non_Active_Receptor_concentration(1); Active_Receptor_concentration(1)];

% Set ODE options for step size

options = odeset('MaxStep', 0.01, 'RelTol', 1e-6, 'AbsTol', 1e-8);

t = [];

y = [];

for i = 1:numel(PhaseTimes)-1

PhaseTime = [PhaseTimes(i) PhaseTimes(i+1)];

% Define functions for forward and backward reaction rates

Kf_L = @(t) calculate_kf(t, PhaseTime, Kf_Max, RLC(i), TauKf_ON, TauKf_OFF);

Kb = @(t) calculate_kb(t, PhaseTime, Kb_Max, Kb_Min, RLC(i), TauKb_ON, TauKb_OFF);

[T,Y] = ode45(@(t,y)ode_LR(t,y,Kf_L,Kb),PhaseTime,initial_conditions, options);

initial_conditions = Y(end,:);

t = [t;T];

y = [y;Y];

end

% Extract the concentrations

Non_Active_Receptor_concentration = y(:, 1);

Active_Receptor_concentration = y(:, 2);

% Calculate forward and backward reaction rates over time

Kf_values = arrayfun(Kf_L, t);

Kb_values = arrayfun(Kb, t);

% Plot active and non-active receptor concentrations

figure;

plot(t, Non_Active_Receptor_concentration, 'r', 'DisplayName', 'Non-Active Receptor Concentration');

hold on;

plot(t, Active_Receptor_concentration, 'b', 'DisplayName', 'Active Receptor Concentration');

legend;

xlabel('Time');

ylabel('Concentration');

title('Receptor Concentrations');

hold off;

% Plot forward and backward reaction rates

figure;

plot(t, Kf_values, 'k', 'DisplayName', 'Forward Reaction Rate');

hold on;

plot(t, Kb_values, 'b', 'DisplayName', 'Backward Reaction Rate');

legend;

xlabel('Time');

ylabel('Reaction Rate');

title('Reaction Rates');

hold off;

% Calculate phase results

PhaseResults = arrayfun(@(i) ...

calculate_phase(t, Active_Receptor_concentration, PhaseTimes(:,i:i+1), RLC(i)) ...

, 1:(size(PhaseTimes, 2)-1), 'UniformOutput', false);

PhaseResults = vertcat(PhaseResults{:}); % Concatenate results

% Calculate maximum forward reaction rate using element-wise division

Kf_LMax = Kf_Max * (RLC ./ (RLC + 1));

% Write Phase results to CSV

for i = 1:size(PhaseResults, 1)

PhaseTable = struct2table(PhaseResults(i));

writetable(PhaseTable, ['Phase', num2str(i), '.csv']);

end

% Nested function for calculating Kf_L based on time

function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)

% Calculate maximum Kf_L based on RLC values using element-wise division

Kf_LMax = Kf_Max * (RLC ./ (1 + RLC));

% Default Kf_L to the maximum value of the first phase

Kf_L = Kf_LMax(1);

% Number of phases

num_phases = numel(RLC);

% Iterate through each phase to determine the current phase based on time

for j = 1:num_phases

T_start = PhaseTimes(j); % Start time of the current phase

if j < num_phases

T_end = PhaseTimes(j + 1); % End time of the current phase

else

T_end = inf; % Handle last phase separately

end

% Check if the current time t is within the current phase

if t >= T_start && t < T_end

if j == 1

% For the first phase, use the maximum value directly

Kf_L = Kf_LMax(j);

else

% Time at the end of the previous phase

prev_end = PhaseTimes(j-1);

% Determine if RLC is increasing or decreasing and calculate Kf_L

if RLC(j) > RLC(j - 1)

% Increasing phase

Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (t - prev_end));

Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));

elseif RLC(j) < RLC(j - 1)

% Decreasing phase

Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (t - prev_end));

Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));

end

return; % Exit the function after finding the correct phase

end

% Nested function for calculating Kb based on time

function Kb = calculate_kb(t, PhaseTimes, Kb_Max, Kb_Min, RLC, TauKb_ON, TauKb_OFF)

% Default to Kb_Max

Kb = Kb_Max;

% If all RLC values are the same, set Kb to Kb_Max and return

if all(RLC == RLC(1))

Kb = Kb_Max;

return;

end

% Iterate through each phase based on time

for j = 1:numel(RLC)

T_start = PhaseTimes(j); % Start time of the current phase

if j < numel(PhaseTimes)

T_end = PhaseTimes(j + 1); % End time of the current phase

else

T_end = inf; % For the last phase

end

% Check if current time is within the current phase

if t >= T_start && t < T_end

if j == 1

% In the first phase, set Kb to Kb_Max

Kb = Kb_Max;

else

prev_end = PhaseTimes(j-1); % Time at the end of the previous phase

% Handle increasing or decreasing RLC values

if RLC(j) > RLC(j - 1)

% Increasing phase

Kb_end = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));

Kb = Kb_Min + (Kb_end - Kb_Min) * exp(TauKb_OFF * (t - prev_end));

elseif RLC(j) < RLC(j - 1)

% Decreasing phase

Kb_end = Kb_Max - (Kb_Max - Kb_Min) * exp(TauKb_ON * (t - prev_end));

Kb = Kb_Min + (Kb_end - Kb_Min) * exp(TauKb_OFF * (t - prev_end));

end

return; % Exit once the correct phase is found

end

% ODE for calculation of the derivative

function dydt = ode_LR(t, y, Kf_L, Kb)

Non_Active_Receptor_concentration = y(1);

Active_Receptor_concentration = y(2);

% Call the function handles for Kf_L and Kb

dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;

dActiveReceptor_dt = Kf_L(t) * Non_Active_Receptor_concentration - Kb(t) * Active_Receptor_concentration;

dydt = [dNonActiveReceptor_dt; dActiveReceptor_dt];

end

% Results of all the Phase Values

function result = calculate_phase(t, Active_Receptor_concentration, PhaseTimes, RLC)

T_start = PhaseTimes(:,1);

T_end = PhaseTimes(:,2);

Phase_indices = (t >= T_start & t <= T_end);

Phase_concentration = Active_Receptor_concentration(Phase_indices);

Phase_time = t(Phase_indices);

% Calculate peak and steady-state values

[Rpeak, Tpeak, peak_type] = findpeak(Phase_time, Phase_concentration);

Rss = mean(Active_Receptor_concentration(t >= (T_end - 10) & t <= T_end));

% Calculate the T50 (time to reach half of peak value)

half_peak_value = (Rss + Rpeak) / 2;

[~, idx_50_percent] = min(abs(Phase_concentration - half_peak_value));

T50 = Phase_time(idx_50_percent) - Tpeak;

% Calculate other metrics

ratio_Rpeak_Rss = Rpeak / Rss;

Delta = Rpeak - Rss;

L = RLC;

% Package results in a struct

result.Rpeak = Rpeak;

result.Rss = Rss;

result.Tpeak = Tpeak;

result.T50 = T50;

result.ratio_Rpeak_Rss = ratio_Rpeak_Rss;

result.Delta = Delta;

result.L = L;

result.peak_type = peak_type;

end

% Customize function for Calculating the Peak values and time

function [Rpeak, Tpeak, peak_type] = findpeak(time, concentration)

% Compute the derivative

dt = diff(time);

dy = diff(concentration);

derivative = dy ./ dt;

% Find zero-crossings of the derivative

zero_crossings = find(diff(sign(derivative)) ~= 0);

% Initialize output variables

Rpeak = NaN;

Tpeak = NaN;

peak_type = 'None';

% Check if there are zero crossings

if ~isempty(zero_crossings)

% Identify peaks and troughs by examining the sign changes

max_indices = zero_crossings(derivative(zero_crossings) > 0 & derivative(zero_crossings + 1) < 0);

min_indices = zero_crossings(derivative(zero_crossings) < 0 & derivative(zero_crossings + 1) > 0);

% Check if there are maxima or minima

if ~isempty(max_indices) || ~isempty(min_indices)

if ~isempty(max_indices) && ~isempty(min_indices)

% Find the highest maximum

[Rpeak_max, maxIdx] = max(concentration(max_indices));

% Find the lowest minimum

[Rpeak_min, minIdx] = min(concentration(min_indices));

% Determine whether the highest maximum is greater or the lowest minimum is smaller

if Rpeak_max >= abs(Rpeak_min)

Rpeak = Rpeak_max;

Tpeak = time(max_indices(maxIdx));

peak_type = 'High';

else

Rpeak = Rpeak_min;

Tpeak = time(min_indices(minIdx));

peak_type = 'Low';

end

% If only maxima exist

elseif ~isempty(max_indices)

[Rpeak, maxIdx] = max(concentration(max_indices));

Tpeak = time(max_indices(maxIdx));

peak_type = 'High';

% If only minima exist

elseif ~isempty(min_indices)

[Rpeak, minIdx] = min(concentration(min_indices));

Tpeak = time(min_indices(minIdx));

peak_type = 'Low';

end

Torsten 2024-9-23，10:25

编辑：Torsten 2024-9-23，10:27

You ask me questions about the code that you made. How should I know better than you do ?

I set the initial condition for the new phase to the end value of the last phase. Maybe that's the problem.

for i = 1:numel(PhaseTimes)-1

PhaseTime = [PhaseTimes(i) PhaseTimes(i+1)];

% Define functions for forward and backward reaction rates

Kf_L = @(t) calculate_kf(t, PhaseTime, Kf_Max, RLC(i), TauKf_ON, TauKf_OFF);

Kb = @(t) calculate_kb(t, PhaseTime, Kb_Max, Kb_Min, RLC(i), TauKb_ON, TauKb_OFF);

[T,Y] = ode45(@(t,y)ode_LR(t,y,Kf_L,Kb),PhaseTime,initial_conditions, options);

initial_conditions = Y(end,:);

t = [t;T];

y = [y;Y];

end

Sam Chak 2024-9-23，11:33

Hi @Ehtisham

If you are unfamiliar with the material, it is advisable not to use AI-generated code (or code inherited from your PhD seniors or kind professor) to describe what you want to achieve, as the code may contain flaws of which you are unaware.

Instead, it will be clearer if you use the mathematical model of the RLC circuit to explain the problem and the parameters you wish to vary throughout the simulation. Click

to insert image of the RLC math model.

Additionally, if the output responses are not as expected, consider applying a "divide-and-conquer" approach by simulating the RLC circuit in phases, individually. For example, if Phase 1 is successful, then introduce Phase 2 into the code, as Phase 2 depends on the results of Phase 1. Simulating all phases at once can complicate troubleshooting due to the various intertwined factors, including coding, human errors, and modeling issues.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Steven Lord 2024-9-20，12:56

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2154185-how-can-i-fix-ode45-returning-a-straight-line-when-using-specific-tspan-points-instead-of-showing#answer_1520105

在 MATLAB Online 中打开

Your ode_LR function is defined to accept four input arguments, the two required arguments plus two additional parameters.

function dydt = ode_LR(t, y, Kf_L, Kb)

How are you using it, based on what's displayed in the error message?

Error in SimfileNPhase>@(t,y)ode_LR(t,y,RLC(i)) (line 20)
  [t, y] = ode45(@(t,y)ode_LR(t,y,RLC(i)),PhaseTime,initial_conditions, options);

Not only are you not calling it with enough input arguments (the anonymous function calls it with three inputs) the comment inside ode_LR implies that the Kf_L and Kb input arguments are function handles.

% Call the function handles for Kf_L and Kb
    dNonActiveReceptor_dt = -Kf_L(t) * Non_Active_Receptor_concentration + Kb(t) * Active_Receptor_concentration;

Is RLC(i), the value that you pass into ode_LR as the third input argument, a function handle? No, it is not. It's a numeric array. And t (the value that ode45 passes into your function as the first required input) is not guaranteed to be a valid index into a numeric array.