How to fix issue with error using surf: Z must be a matrix, not a scalar or vector

How can I fix the code to avoid the error: error using surf: Z must be a matrix, not a scalar or vector

I tried converting it to matrix format using cell2mat, but then I can't refer to the correct position of terms to plot on the surface. Is there an alternative way?

numODEs =3;
numParams=3;

% Parameters
parameters = rand(1, numParams);  % Initialize parameters with random values

% Initial conditions
initial_conditions = rand(1, numODEs);

% Time span for simulation
tspan = [0 10];

% Set variation range based on a percentage of original parameters
min_variation_factor = 0.1;
max_variation_factor = 10;
variation_range = linspace(min_variation_factor, max_variation_factor, 5);

% Store results
balance_results = cell(numParams, numParams);

% Explore parameter balance for all combinations
for i = 1:numParams
    for j = 1:numParams
        % Perturb parameters based on variation_range
        perturbed_params = parameters;
        perturbed_params(i) = perturbed_params(i) * variation_range(1);
        perturbed_params(j) = perturbed_params(j) * variation_range(2);

        % Solve ODE system with perturbed parameters using ode15s
        [~, Y_perturbed] = ode15s(@ode_system, tspan, initial_conditions, [], perturbed_params);

        % Check if the solution compensates for variations
        balance_results{i, j} = max(Y_perturbed, [], 1); % Adjust based on your specific measure
        
    end
end

balance_results{1, 3}(:, 1)

ans = 0.8496

% Visualize parameter balance
    figure;
    for k = 1:numODEs
        subplot(ceil(sqrt(numODEs)), ceil(sqrt(numODEs)), k);
        colormap(parula(numParams^2));
        for i = 1:numParams
           for j = 1:numParams
               surf(variation_range(1), variation_range(2), balance_results{i, j}(:, k)', 'FaceColor', 'interp', 'EdgeColor', 'none');
               hold on;
           end
       end
       xlabel('Variation in Parameter 1');
       ylabel('Variation in Parameter 2');
        zlabel(['Max Output for ODE ', num2str(k)]);
        title(['Balance Exploration for ODE ', num2str(k)]);
        view(-45, 30); % Adjust the view for better visibility
        hold off;
    end

Error using surf
Z must be a matrix, not a scalar or vector.

function dydt = ode_system(t, y, parameters)
    % ODE system representing binding reactions (modify based on your specific system)
    dydt = zeros(size(y));
    
% Parameters
    k1 = parameters(1);
    k2 = parameters(2);
    k3 = parameters(3);

    dydt(1) = -k1 * y(1) * y(2);  % Example binding reaction
    dydt(2) = k1 * y(1) * y(2) - k2 * y(2);  % Example binding and unbinding reaction
    dydt(3) = k2 * y(2) - k3 * y(3);  % Example unbinding reaction
    % Define the ODEs (modify based on your specific system)
    % Example: dydt(1) = -parameters(1) * y(1) * y(2);

    % Modify the code based on the structure of your ODE system
end

[Dyuman Joshi]

How do you plan to plot a surface where the data you have is scalar? You can't even plot a line with a single point.

What are you trying to do? What is the objective here?

[Rebeca Hannah Oliveira]

Thank you! My goal is to get something like in the code below (which was specific for three ODEs and 3 parameters), but generalized to any number n of ODEs and parameters p. I want to generate n figures, each with p subplots that are surface plots showing what happens to the output of the ODEs when two parameters are varied within a range.

% function explore_parameter_balance_binding_reactions()
    % Number of ODEs and parameters
    numODEs = 3;
    numParams = 3;

    % Parameters
    parameters = rand(1, numParams);  % Initialize parameters with random values

    % Initial conditions
    initial_conditions = rand(1, numODEs);

    % Time span for simulation
    tspan = [0 10];

    % Set variation range based on a percentage of original parameters
    min_variation_factor = 0.1;
    max_variation_factor = 10;
    variation_range = linspace(min_variation_factor, max_variation_factor, 5);

    % Store results
    balance_results_1_2 = zeros(length(variation_range), length(variation_range), numODEs);
    balance_results_1_3 = zeros(length(variation_range), length(variation_range), numODEs);
    balance_results_2_3 = zeros(length(variation_range), length(variation_range), numODEs);

    % Explore parameter balance for all combinations
    for i = 1:length(variation_range)
        for j = 1:length(variation_range)
            % Perturb parameters based on variation_range
            parameters_perturbed_1_2 = parameters .* [variation_range(i), variation_range(j), 1];
            parameters_perturbed_1_3 = parameters .* [variation_range(i), 1, variation_range(j)];
            parameters_perturbed_2_3 = parameters .* [1, variation_range(i), variation_range(j)];

            % Solve ODE system with perturbed parameters using ode15s
            [~, Y_perturbed_1_2] = ode15s(@ode_system, tspan, initial_conditions, [], parameters_perturbed_1_2);
            [~, Y_perturbed_1_3] = ode15s(@ode_system, tspan, initial_conditions, [], parameters_perturbed_1_3);
            [~, Y_perturbed_2_3] = ode15s(@ode_system, tspan, initial_conditions, [], parameters_perturbed_2_3);

            % Check if the solution compensates for variations
            balance_results_1_2(i, j, :) = max(Y_perturbed_1_2, [], 1); % Adjust based on your specific measure
            balance_results_1_3(i, j, :) = max(Y_perturbed_1_3, [], 1);
            balance_results_2_3(i, j, :) = max(Y_perturbed_2_3, [], 1);
        end
    end

    % Visualize parameter balance for parameters 1 and 2
    figure;
    for k = 1:numODEs
        subplot(ceil(sqrt(numODEs)), ceil(sqrt(numODEs)), k);
        surf(variation_range, variation_range, squeeze(balance_results_1_2(:, :, k))');
        xlabel('Variation in Parameter 1');
        ylabel('Variation in Parameter 2');
        zlabel(['Max Output for ODE ', num2str(k)]);
        title(['Balance Exploration for ODE ', num2str(k), ' (Params 1 and 2)']);
    end
    sgtitle('Parameter Balance Exploration for Binding Reactions (Params 1 and 2)');

    % Visualize parameter balance for parameters 1 and 3
    figure;
    for k = 1:numODEs
        subplot(ceil(sqrt(numODEs)), ceil(sqrt(numODEs)), k);
        surf(variation_range, variation_range, squeeze(balance_results_1_3(:, :, k))');
        xlabel('Variation in Parameter 1');
        ylabel('Variation in Parameter 3');
        zlabel(['Max Output for ODE ', num2str(k)]);
        title(['Balance Exploration for ODE ', num2str(k), ' (Params 1 and 3)']);
    end
    sgtitle('Parameter Balance Exploration for Binding Reactions (Params 1 and 3)');

    % Visualize parameter balance for parameters 2 and 3
    figure;
    for k = 1:numODEs
        subplot(ceil(sqrt(numODEs)), ceil(sqrt(numODEs)), k);
        surf(variation_range, variation_range, squeeze(balance_results_2_3(:, :, k))');
        xlabel('Variation in Parameter 2');
        ylabel('Variation in Parameter 3');
        zlabel(['Max Output for ODE ', num2str(k)]);
        title(['Balance Exploration for ODE ', num2str(k), ' (Params 2 and 3)']);
    end
    sgtitle('Parameter Balance Exploration for Binding Reactions (Params 2 and 3)');

% end

function dydt = ode_system(t, y, parameters)
    % ODE system representing binding reactions (modify based on your specific system)
    dydt = zeros(size(y));

    % Parameters
    k1 = parameters(1);
    k2 = parameters(2);
    k3 = parameters(3);

    % Define the ODEs (modify based on your specific system)
    dydt(1) = -k1 * y(1) * y(2);  % Example binding reaction
    dydt(2) = k1 * y(1) * y(2) - k2 * y(2);  % Example binding and unbinding reaction
    dydt(3) = k2 * y(2) - k3 * y(3);  % Example unbinding reaction
end

Answer 1

[Pratyush]

Hi Rebecca,

I understand that you are getting an error using the "surf" function in your script.

You could use "meshgrid" function to resolve the error. To fix this, you need to create a meshgrid of "variation_range(1)" and "variation_range(2)" to use as the X and Y inputs for the "surf" function. Then you can use the "meshgrid" function to create the grid and then plot the surfaces accordingly. Here's how you can modify your code:

% Visualize parameter balance
figure;
for k = 1:numODEs
    subplot(ceil(sqrt(numODEs)), ceil(sqrt(numODEs)), k);
    colormap(parula(numParams^2));
    [X, Y] = meshgrid(variation_range, variation_range);
    for i = 1:numParams
        for j = 1:numParams
            Z = balance_results{i, j}(:, k)';
            surf(X, Y, reshape(Z, size(X)), 'FaceColor', 'interp', 'EdgeColor', 'none');
            hold on;
        end
    end
    xlabel('Variation in Parameter 1');
    ylabel('Variation in Parameter 2');
    zlabel(['Max Output for ODE ', num2str(k)]);
    title(['Balance Exploration for ODE ', num2str(k)]);
    view(-45, 30); % Adjust the view for better visibility
    hold off;
end

I hope this should resolve the error.

[Dyuman Joshi]

There's an error when I try to run the code with your suggestion -

numODEs =3;
numParams=3;

% Parameters
parameters = rand(1, numParams);  % Initialize parameters with random values

% Initial conditions
initial_conditions = rand(1, numODEs);

% Time span for simulation
tspan = [0 10];

% Set variation range based on a percentage of original parameters
min_variation_factor = 0.1;
max_variation_factor = 10;
variation_range = linspace(min_variation_factor, max_variation_factor, 5);

% Store results
balance_results = cell(numParams, numParams);

% Explore parameter balance for all combinations
for i = 1:numParams
    for j = 1:numParams
        % Perturb parameters based on variation_range
        perturbed_params = parameters;
        perturbed_params(i) = perturbed_params(i) * variation_range(1);
        perturbed_params(j) = perturbed_params(j) * variation_range(2);

        % Solve ODE system with perturbed parameters using ode15s
        [~, Y_perturbed] = ode15s(@ode_system, tspan, initial_conditions, [], perturbed_params);

        % Check if the solution compensates for variations
        balance_results{i, j} = max(Y_perturbed, [], 1); % Adjust based on your specific measure
        
    end
end


balance_results{1, 3}(:, 1)

ans = 0.3748

% Visualize parameter balance
figure;
for k = 1:numODEs
    subplot(ceil(sqrt(numODEs)), ceil(sqrt(numODEs)), k);
    colormap(parula(numParams^2));
    [X, Y] = meshgrid(variation_range, variation_range);
    for i = 1:numParams
        for j = 1:numParams
            Z = balance_results{i, j}(:, k)';
            surf(X, Y, reshape(Z, size(X)), 'FaceColor', 'interp', 'EdgeColor', 'none');
            hold on;
        end
    end
    xlabel('Variation in Parameter 1');
    ylabel('Variation in Parameter 2');
    zlabel(['Max Output for ODE ', num2str(k)]);
    title(['Balance Exploration for ODE ', num2str(k)]);
    view(-45, 30); % Adjust the view for better visibility
    hold off;
end

Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.

function dydt = ode_system(t, y, parameters)
    % ODE system representing binding reactions (modify based on your specific system)
    dydt = zeros(size(y));
    
% Parameters
    k1 = parameters(1);
    k2 = parameters(2);
    k3 = parameters(3);

    dydt(1) = -k1 * y(1) * y(2);  % Example binding reaction
    dydt(2) = k1 * y(1) * y(2) - k2 * y(2);  % Example binding and unbinding reaction
    dydt(3) = k2 * y(2) - k3 * y(3);  % Example unbinding reaction
    % Define the ODEs (modify based on your specific system)
    % Example: dydt(1) = -parameters(1) * y(1) * y(2);

    % Modify the code based on the structure of your ODE system
end

[Rebeca Hannah Oliveira]

Pratyush, thank you. It is returning an error with reshape, though. Can you help again?