File:
Transclude of bifurcations.m(It will require some matlab packs, like “symbolik toolbox”)
Instructions:
- Copy the script.
- Paste it into MATLAB (It does work for sure in MATLAB Online, as of writing
2024/10/01) - Run it and see what it does.
Change the my_system function and insert the sytem you need:
- The
varsarray contains the variables of the system. - The
subs_matrix contains the fix parameters of the system. - The
paramis the parameter of the system. - The
tspantime frame used for the simulation. - The
initial_conditionsare the initial_conditions. - The
param_valueis the value of the parameter used for the simulation.
Also I have defined many systems to help you understand how do properly define a new one, so that the program can succefully run, ==you can find them in the SYSTEMS section==.
Script:
% MATLAB sucks: TOO many warnings
% MATLAB mixed feelings: however the way to suppress them is explicit
%#ok<*DEFNU>
%#ok<*NOPTS>
%#ok<*NOPRT>
%#ok<*ISCL>
close all
clear
clc
MAIN -------------------------------------------------------------------
[system, vars, subs_, param, param_range, ...
tspan, initial_conditions, param_value] = data();
system
steady_states = solve(rhs(system) == 0, vars);
steady_states = reshape_solutions(steady_states);
fprintf('num_of_steady_states: %i\n\n', length(steady_states))
for i=1:length(steady_states)
fprintf('steady_state{%i}:\n', i);
disp(steady_states{i})
end
J = jacobian(rhs(system), vars)
if ~isempty(subs_)
sys = subs(rhs(system), subs_(:,1), subs_(:,2));
J = subs(J, subs_(:,1), subs_(:,2))
else
sys = rhs(system)
end
analyze_each_steady_state(J, steady_states, vars, subs_, param, param_range)
function analyze_each_steady_state(J, steady_states, vars, subs_, param, param_range)
disp('--- FUNCTION: analyze_each_steady_state ---')
for i=1:length(steady_states)
if ~isempty(subs_)
ss = subs(steady_states{i}, subs_(:,1), subs_(:,2))
else
ss = steady_states{i}
end
J_lin = subs(J, vars, ss)
bifurcation_values = find_bifurcations_values(J_lin, param, param_range)
check_stability(J_lin, ss, bifurcation_values, param)
disp('-----------------------------------------')
end
end
tspan
initial_conditions
param_value
% syms t
sys_ = subs(sys, param, param_value)
ht = @(t, xyz) double(subs(sys_, vars, xyz));
[t, xyz] = ode_func(ht, tspan, initial_conditions, ode_options);
figure
hold on
if length(vars) == 1
x = xyz(:,1);
plot(t, x)
plot(t(1), x(1), 'go','MarkerFaceColor','g')
plot(t(end), x(end), 'bo','MarkerFaceColor','b')
xlabel('t')
ylabel('x(t)')
legend('x(t)', 'START', 'END')
% disp("TODO")
elseif length(vars) == 2
x = xyz(:,1);
y = xyz(:,2);
plot(x, y)
plot(x(1), y(1), 'go','MarkerFaceColor','g')
plot(x(end), y(end), 'bo','MarkerFaceColor','b')
dx = diff(x)
dy = diff(y)
quiver(x(1:end-1), y(1:end-1), dx, dy, 0)
xlabel('x(t)')
ylabel('y(t)')
legend('', 'START', 'END', '(x, y)')
elseif length(vars) == 3
x = xyz(:,1);
y = xyz(:,2);
z = xyz(:,3);
plot3(x, y, z)
plot3(x(1), y(1), z(1), 'go','MarkerFaceColor','g')
plot3(x(end), y(end), z(end), 'bo','MarkerFaceColor','b')
xlabel('x(t)')
ylabel('y(t)')
legend('(x, y, z)', 'START', 'END')
elseif length(vars) >= 4
disp("TODO")
% to try
legend_array = string(zeros(length(vars),1));
for i=1:length(vars)
x_i = xyz(:,i);
plot(t, x_i);
legend_array(i) = strcat('x', string(i));
end
xlabel('t')
ylabel('x_i(t)')
legend(legend_array)
end
hold off
FUNCTIONS --------------------------------------------------------------
function [out] = reshape_solutions(steady_states)
parameters = fieldnames(steady_states);
if isempty(parameters)
out = num2cell(steady_states);
return
end
first_param_size = size(steady_states.(parameters{1}));
assert(first_param_size(2) == 1)
for i = 1:length(parameters)
assert(all(size(steady_states.(parameters{i}))==first_param_size))
end
num_of_sol = size(steady_states.x, 1);
out = cell(num_of_sol,1);
for i=1:num_of_sol
temp = [];
for j=1:length(parameters)
temp = [temp; steady_states.(parameters{j})(i)];
end
out{i} = temp;
end
end
function [bifurcation_values] = find_bifurcations_values(J_lin, param, param_range)
eigenvalues = eig(J_lin);
bifurcation_values = double([]);
for i=1:length(eigenvalues)
sym_vars = symvar(real(eigenvalues(i)));
num_sym_vars = length(sym_vars);
if num_sym_vars ~= 1
continue
end
solutions = vpasolve(real(eigenvalues(i))==0, param, param_range);
for j=1:length(solutions)
new_value = double(solutions(j));
if ismember(new_value, bifurcation_values)
continue
end
bifurcation_values(end+1) = new_value;
end
end
end
function check_stability(J_lin, ss, bifurcation_values, param)
param_name = string(symvar(param));
for i=1:length(bifurcation_values)
if bifurcation_values(i) == 0
k_plus = 0.0001;
k_minus = -0.0001;
else
k_plus = bifurcation_values(i)*1.0001;
k_minus = bifurcation_values(i)*0.9999;
end
J_rp = subs(J_lin, param, k_plus);
J_rm = subs(J_lin, param, k_minus);
eigs_kp = double(eig(J_rp));
eigs_km = double(eig(J_rm));
num_of_negative_eigs = sum(real(eigs_km) < 0);
num_of_positive_eigs = sum(real(eigs_km) > 0);
num_of_null_eigs = sum(real(eigs_km) == 0);
assert(num_of_null_eigs == 0)
is_imaginary = any(imag(subs(ss, param, k_minus)) ~= 0);
str_msg = "";
if is_imaginary
str_msg = "(IMAGINARY steady state)";
end
assert(num_of_null_eigs == 0)
fprintf('For %s < %.3f : %s \n\t%i negative eigenvalues, %i positive eigenvalues\n', ...
param_name, ...
bifurcation_values(i), ...
str_msg, ...
num_of_negative_eigs, ...
num_of_positive_eigs ...
);
num_of_negative_eigs = sum(real(eigs_kp) < 0);
num_of_positive_eigs = sum(real(eigs_kp) > 0);
num_of_null_eigs = sum(real(eigs_kp) == 0);
assert(num_of_null_eigs == 0)
is_imaginary = any(imag(subs(ss, param, k_plus)) ~= 0);
str_msg = "";
if is_imaginary
str_msg = "(IMAGINARY steady state)";
end
assert(num_of_null_eigs == 0)
fprintf('For %s > %.3f : %s \n\t%i negative eigenvalues, %i positive eigenvalues\n', ...
param_name, ...
bifurcation_values(i), ...
str_msg, ...
num_of_negative_eigs, ...
num_of_positive_eigs ...
);
end
end