求助：Matlab分岔图代码运行耗时久且无结果问题排查

Here are the critical issues in your code and targeted fixes to resolve the long runtime and lack of output:

Key Problems in Original Code #

• Inefficient array appending: Using results = [results; ...] copies the entire array every iteration, leading to exponential slowdown as the dataset grows.
• Storing unnecessary time points: Bifurcation diagrams only require points from the system's attractor (steady state or periodic orbit), not every time step from t=0 to tmax.
• Redundant ODE definition: Defining the ODE system inside the loop wastes computation time re-creating the function each iteration.

Fixed MATLAB Code #

clc; close all; clear;

% System parameters
r1 = 0.36; a = 0.1; k1 = 0.36; k2 = 0.48; r2 = 1; deltav = 0.2; deltaz = 0.036;

% Simulation settings
t_transient = 200; % Time to let transients decay
t_attractor = 100; % Time to collect attractor points
b_values = 0:0.1:30;
num_b = length(b_values);

% Preallocate results using cell array to avoid size issues
results_cell = cell(num_b, 1);

% Initial conditions
x0 = 0.5; y0 = 0; v0 = 0.3; z0 = 0.1; 

% Define ODE system once outside the loop (pass b as a parameter)
ode_system = @(t, y, b) [
    r1 * y(1) * (1 - y(1)) - a * y(1) * y(3) - k1 * y(1) * y(4);
    a * y(1) * y(3) - k2 * y(2) * y(4) - y(2);
    b * y(2) - a * y(1) * y(3) - deltav * y(3);
    r2 * y(2) * y(4) - deltaz * y(4)
];

% Solver options (balance speed and accuracy)
opts = odeset('RelTol', 1e-6, 'AbsTol', 1e-8);

for idx = 1:num_b
    b = b_values(idx);
    
    % First run to eliminate transients
    [~, sol_transient] = ode45(@(t,y) ode_system(t,y,b), [0, t_transient], [x0, y0, v0, z0], opts);
    final_transient = sol_transient(end, :);
    
    % Collect points from the attractor
    [t_attract, sol_attract] = ode45(@(t,y) ode_system(t,y,b), [0, t_attractor], final_transient, opts);
    
    % Store b values and corresponding x values
    results_cell{idx} = [b * ones(length(t_attract), 1), sol_attract(:, 1)];
    
    % Update initial conditions for next b value
    x0 = sol_attract(end, 1);
    y0 = sol_attract(end, 2);
    v0 = sol_attract(end, 3);
    z0 = sol_attract(end, 4);
end

% Combine cell array into a single matrix for plotting
results = vertcat(results_cell{:});

% Generate bifurcation diagram
figure;
plot(results(:, 1), results(:, 2), '.', 'MarkerSize', 1);
xlabel('b');
ylabel('x');
title('Bifurcation Diagram');
grid on;

Additional Optimizations #

• Steady-state only: If your system settles to a fixed point, replace the attractor collection step with storing just the final value of sol_attract (one point per b instead of multiple time steps).
• Adjust solver tolerance: Loosen RelTol and AbsTol (e.g., to 1e-4) for faster computation if minor accuracy tradeoffs are acceptable.
• Reduce t_attractor: Shorten this value if your system reaches its attractor quickly.

内容的提问来源于stack exchange，提问作者ZSN