When i remove the element wise divison of the RLC beacuse RLC is single array it gives me this error Index exceeds the number of array elements. Index must not exceed 1.

Question

Ehtisham 2024-9-4，9:16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2149989-when-i-remove-the-element-wise-divison-of-the-rlc-beacuse-rlc-is-single-array-it-gives-me-this-error

评论： Sam Chak 2024-9-5，13:07

Need your guidance

% Define parameters
Kf_Max = 100;  % maximum forward rate
RLC = [0.1, 0.5, 10, 1];  % RLC values
TauKf_ON = -0.1; % TauKf_ON
TauKf_OFF = -0.1;  % TauKf_OFF
Kb_Max = 50;  % maximum backward rate  
PhaseTimes = [0, 100, 200, 300, 400];% phase times (each row defines a phase)
timespan = [0,  400];  % timespan for the simulation
% 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,  PhaseTimes, timespan, 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,  PhaseTimes, timespan, Non_Active_Receptor_concentration, Active_Receptor_concentration)
    % Define functions for forward and backward reaction rates
    Kf_L = @(t) calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF);
    Kb = @(t) calculate_kb( Kb_Max);
    % 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-4, 'AbsTol', 1e-6);
    % Solve the ODE system
    [t, y] = ode45(@(t, y) ode_LR(t, y, Kf_L, Kb), timespan, initial_conditions, options);
    % 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, 'c', '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
    Kf_LMax = Kf_Max * (RLC / (1+ RLC));
    % Write Phase results to CSV
    for i = 1:size(PhaseResults, 1)
        PhaseTable = struct2table(PhaseResults(i));
        writetable(PhaseTable, ['Phase', num2str(i), '.csv']);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kf_L = calculate_kf(t, PhaseTimes, Kf_Max, RLC, TauKf_ON, TauKf_OFF)
    % Calculate the maximum Kf_L for each phase based on RLC values
    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 based on the length of RLC
    num_phases = numel(RLC);
    % Iterate through each phase to determine the current phase based on time t
    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);
                % For phases other than the last one
                if j < num_phases
                    % Peak condition: RLC value is higher in the current phase compared to its neighbors
                    if RLC(j - 1) < RLC(j) && RLC(j) > RLC(j + 1)
                        Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (T_end - T_start));
                        Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                    % Increasing RLC condition: 
                    elseif RLC(j - 1) < RLC(j) && RLC(j) <= RLC(j + 1)
                        Kf_L = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (T_end - T_start));
                       
                    % Decreasing RLC condition: RLC value is lower in the current phase compared to neighbors
                    elseif RLC(j - 1) > RLC(j) && RLC(j) < RLC(j + 1)
                        Kf_L = Kf_LMax(j) + (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                    % Flat or decreasing RLC condition
                    elseif RLC(j - 1) > RLC(j) && RLC(j) >= RLC(j + 1)
                        Kf_L = Kf_LMax(j) + (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                    end
                else
                    % For the last phase, consider only the previous phase's value and time
                    if RLC(j - 1) < RLC(j)
                        Kf_end = Kf_LMax(j) - (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (T_end - T_start));
                        Kf_L = Kf_LMax(j) + (Kf_end - Kf_LMax(j - 1)) * exp(TauKf_OFF * (t - prev_end));
                    elseif RLC(j - 1) > RLC(j)
                        Kf_L = Kf_LMax(j) + (Kf_LMax(j) - Kf_LMax(j - 1)) * exp(TauKf_ON * (t - prev_end));
                    end
                end
            end
            return;
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Kb = calculate_kb( Kb_Max)
    Kb = Kb_Max;  % Default to the maximum value
end
% function Kb = calculate_kb(t, PhaseTimes, Kb_Max, RLC)
%     % Default value for Kb is set to Kb_Max
%     Kb = Kb_Max; 
% 
%     % Check if all RLC values are the same
%     if all(RLC == RLC(1))
%         % If all RLC values are the same, set Kb to Kb_Max
%         Kb = Kb_Max; 
%         return;
%     end
% 
%     % Loop through each phase interval
%     for j = 1:numel(RLC)
%         % Define the start and end times for the current phase
%         T_start = PhaseTimes(j);
%         T_end = PhaseTimes(j + 1);
% 
%         % Check if current time falls within the current phase interval
%         if t >= T_start && t < T_end
%             % If it's the first phase, Kb is always Kb_Max
%             if j == 1
%                 Kb = Kb_Max;
%             else
%                 % Check if we're not at the last phase
%                 if j < numel(RLC)
%                     % Determine Kb based on the relative values of RLC
%                     % Compare current RLC value with previous and next values
%                     if RLC(j-1) < RLC(j) && RLC(j) > RLC(j+1)
%                         % Local maximum
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) < RLC(j) && RLC(j) <= RLC(j+1)
%                         % Local peak or plateau
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) > RLC(j) && RLC(j) < RLC(j+1)
%                         % Local minimum
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) > RLC(j) && RLC(j) >= RLC(j+1)
%                         % Local dip or plateau
%                         Kb = Kb_Max;
%                     end
%                 else
%                     % Handle the case when j is the last phase
%                     if RLC(j-1) < RLC(j)
%                         % If last value is higher than previous
%                         Kb = Kb_Max;
%                     elseif RLC(j-1) > RLC(j)
%                         % If last value is lower than previous
%                         Kb = Kb_Max;
%                     end
%                 end
%             end
%             return; % Exit function once the appropriate Kb value is set
%         end
%     end
% end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dydt = ode_LR(t, y, Kf_L, Kb)
    Non_Active_Receptor_concentration = y(1);
    Active_Receptor_concentration = y(2);
    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;
% print the Active and NoN Active Receptor concentration 
fprintf('t = %f, Non_Active = %f, Active = %f, dNonActive/dt = %f, dActive/dt = %f\n', t, Non_Active_Receptor_concentration, Active_Receptor_concentration, dNonActiveReceptor_dt, dActiveReceptor_dt);
    dydt = [dNonActiveReceptor_dt; dActiveReceptor_dt];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Rpeak, Tpeak, peak_type] = findpeak(Phase_time, Phase_concentration)
    % Compute the derivative 
    dt = diff(Phase_time);
    dy = diff(Phase_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(Phase_concentration(max_indices));
                % Find the lowest minimum
                [Rpeak_min, minIdx] = min(Phase_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 = Phase_time(max_indices(maxIdx));
                    peak_type = 'High';
                else
                    Rpeak = Rpeak_min;
                    Tpeak = Phase_time(min_indices(minIdx));
                    peak_type = 'Low';
                end
            % If only maxima exist
            elseif ~isempty(max_indices)
                [Rpeak, maxIdx] = max(Phase_concentration(max_indices));
                Tpeak = Phase_time(max_indices(maxIdx));
                peak_type = 'High';
            % If only minima exist
            elseif ~isempty(min_indices)
                [Rpeak, minIdx] = min(Phase_concentration(min_indices));
                Tpeak = Phase_time(min_indices(minIdx));
                peak_type = 'Low';
            end
        end
    end
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Torsten 2024-9-4，11:00

编辑：Torsten 2024-9-4，11:09

在 MATLAB Online 中打开

You must keep the elementwise division. So

Change

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

to

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

What else do you want Kf_LMax(i) to become than Kf_Max * RLC(i)/(1+RLC(i)) ?

And you should use ode15s instead of ode45 - it's way faster.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2024-9-4，13:01

编辑：Torsten 2024-9-4，13:09

在 MATLAB Online 中打开

Do you know what you get from

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

if RLC is a vector ? You get nonsense (at least for your application).

RLC = [1 2 3 4 5];
Kf_Max = 1;
Kf_Max * (RLC / (1 + RLC))
ans = 0.7778
Kf_Max * sum(RLC.*(1+RLC))/sum((1+RLC).^2)
ans = 0.7778

Sam Chak 2024-9-5，13:07

在 MATLAB Online 中打开

Hi @Ehtisham, Essentially, you are stating that the RLC value should "discontinuously jump" to a new value at certain phase times, therefore, the Kf_L value should respond accordingly. However, the original RLC vector contains fewer elements than the PhaseTimes vector. Is this intended?

By the way, did you 'accidentally close' your thread?

https://www.mathworks.com/matlabcentral/answers/2150359-why-the-function-file-is-not-working-on-the-phasetime-time-point-that-define-in-the-script-file?s_tid=srchtitle

%% Modified

RLC_repair = [0, 0.1, 0.5, 10, 1]; % arbitrarily added 0 at the beginning

PhaseTimes = [0, 100, 200, 300, 400];

plot(PhaseTimes, RLC_repair, 'o--'), grid on, grid minor

xlabel('Phase Time'), ylabel('RLC')

%% Original

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

PhaseTimes = [0, 100, 200, 300, 400];

plot(PhaseTimes, RLC, 'o--'), grid on, grid minor

Error using plot
Specify the coordinates as vectors or matrices of the same size, or as a vector and a matrix that share the same length in at least one dimension.

请先登录，再进行评论。

Answer 2

Rahul 2024-9-5，5:43

编辑：Rahul 2024-9-5，6:28

在 MATLAB Online 中打开

Hi Ehtisham,

I understand that you’re trying to solve a series of equations for a given RLC circuit and a set of initial conditions, where an error stating “Index exceeds the number of array elements. Index must not exceed 1” is thrown on program execution.

Specifically, in the function definition of ‘calculate_kft’, since variable ‘Kf_Max’ and division result involving vector ‘RLC’ are scalars themselves, the variable ‘Kf_LMax’ is being assigned a scalar value, which leads to the out of range index error statement:

RLC = [0.1, 0.5, 10, 1];
Kf_Max = 100;
x = (RLC / (1 + RLC))
x = 0.8786
Kf_LMax = Kf_Max * (RLC / (1 + RLC))
Kf_LMax = 87.8561

Further, variable Kb has been assigned a scalar value of ‘Kb_Max’ = 50, inside the function ‘calculate_kb’ definition. Later, these two variables have been assumed to be non-scalars, in the following functions:

calculate_kf: variable ‘num_phases’ stores the size of ‘RLC’ vector, as 3. Hence, while inside the for-loop ‘Iterate through each phase to determine the current phase based on time ‘t’ using an index ‘j’, calculating ‘Kf_end’ and ‘Kf_L’ for ‘phases other than the last one’, ‘Kf_LMax’ is being indexed as a vector, causing the stated error.
ode_LR: variable ‘Kb’ has been indexed using variable ‘t’, whereas earlier ‘Kb’ has been assigned a scalar value:

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;

You can either include element wise division while calculating Kf_LMax as shown below, or make appropriate modifications to the above-mentioned statements.

RLC = [0.1, 0.5, 10, 1];
Kf_Max = 100;
x = (RLC ./ (1 + RLC))
x = 1x4
    0.0909    0.3333    0.9091    0.5000
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>