Mismatched Stability check.
Mostra commenti meno recenti
I have system of two equations. I solved the model - variable k and t for each value of parameter m, between 0 and 1, and checked for stability of solution. I got all solution as stable
% Define the equations as anonymous functions
eq1 = @(vars, m) (sqrt(vars(1))/sqrt(1-vars(1))) + ((vars(2)*(2*(m^3) + vars(2)^3 - 3*vars(2)*((m+1)^2))^2)/(6*vars(2) - 3*((m-vars(2))^2)) + ((m-vars(2))^2)*(m+2*vars(2))*(m/(3*vars(2)) + 2/3)) / ((m^3 - (vars(2)^2)*(3*m+3) + 2*(vars(2)^3))*(m/(3*vars(2)) - (2*m+vars(2))/(3*(2*vars(2)-((m-vars(2))^2))) + 2/3) - vars(2)*((m-vars(2))^2)*(2*m+vars(2))*((2*m/3) + vars(2)/3));
eq2 = @(vars, m) vars(2) - ((sqrt(vars(1))/sqrt(1-vars(1)))*(((m/vars(2))+2)/3) - (1/3)*(2 + (m/vars(2)) + ((2*m+vars(2))/(((m-vars(2))^2)-2*vars(2))))) / ((sqrt(vars(1))/sqrt(1-vars(1)))*(2*(m^3) - 3*((1+m)^2)*vars(2) + vars(2)^3)/(3*(((m-vars(2))^2)-2*vars(2))) - ((2*m+vars(2))/3));
% Range of m values
m_values = linspace(0, 1, 100); % 100 values of m between 0 and 1
% Preallocate arrays to store solutions
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));
stabilities = cell(size(m_values));
% Loop over m values and solve the system of equations for each m
for i = 1:length(m_values)
m = m_values(i);
% System of equations for current m
system_of_equations = @(vars) [eq1(vars, m); eq2(vars, m)];
% Initial guess for [k, t]
initial_guess = [0.75, 1.5];
% Solve the system of equations using fsolve
options = optimoptions('fsolve', 'Display', 'off');
solution = fsolve(system_of_equations, initial_guess, options);
% Store solutions if they meet the criteria
k_solution = solution(1);
t_solution = solution(2);
if k_solution > 0.5 && k_solution < 1 && t_solution > 0
k_solutions(i) = k_solution;
t_solutions(i) = t_solution;
% Compute the Jacobian matrix
J = [diff(eq1([k_solution, t_solution], m), 1, 1), diff(eq1([k_solution, t_solution], m), 1, 2);
diff(eq2([k_solution, t_solution], m), 1, 1), diff(eq2([k_solution, t_solution], m), 1, 2)]
% Calculate eigenvalues of the Jacobian matrix
eigenvalues = eig(J);
% Analyze stability
if all(real(eigenvalues) < 0)
stabilities{i} = 'Stable';
elseif all(real(eigenvalues) == 0)
% Use center manifold method for stability analysis
% Modify this part with your center manifold method implementation
stabilities{i} = 'Center manifold method';
elseif any(real(eigenvalues) > 0)
stabilities{i} = 'Unstable';
end
else
k_solutions(i) = NaN;
t_solutions(i) = NaN;
stabilities{i} = 'NaN';
end
end
% Display solutions and stabilities
disp('Solutions for each m:');
for i = 1:length(m_values)
fprintf('m = %.2f, k = %.4f, t = %.4f, Stability = %s\n', m_values(i), k_solutions(i), t_solutions(i), stabilities{i});
end
Now, the problem arose when I took k as an exogenous variable. When I input the value of k say 0.6945 (In above code, I got one of the solution as for m=0.38 I had k=0.6945 and t =1.82 and it was stable), I got value of t as 1.82 against m=0.38, as should be. But now this solution is unstable. Why? Can anyone help me why I am getting it stable in above code. While it is unstable in below. Please help
% Define the range of m values
m_values = linspace(0, 1, 100); % Adjust the number of points for higher accuracy
% Define the value of k
k = 0.6945;
% Initialize cell arrays to store fixed points and their stability
fixed_points = cell(length(m_values), 1);
stability_info = cell(length(m_values), 1);
% Define the function f(t, m)
f = @(t, m) real(((((k).^(1/2)).*(((m./t)+2)/3) - ((1-k).^(1/2)).*(1/3)*(2 + (m./t) + ((2*m + t)./(((m - t)^2) - 2*t)))))./((((k).^(1/2)).*(2*m^3 - 3*((1+m)^2)*t + t^3)/(3*(((m - t)^2) - 2*t)) - ((1-k).^(1/2)).*((2*m + t)/3))) - t);
% Define initial guesses for fsolve
initial_guesses = [0.1, 0.5, 1.7];
% Loop over the range of m values
for i = 1:length(m_values)
m_val = m_values(i);
% Initialize temporary storage for current m
current_fixed_points = [];
current_stability_info = [];
% Find fixed points using fsolve with multiple initial guesses
for guess = initial_guesses
options = optimset('Display', 'off', 'TolFun', 1e-10, 'TolX', 1e-10);
try
[t_val, fval, exitflag] = fsolve(@(t) f(t, m_val), guess, options);
% Check if the solution is valid and positive
if exitflag > 0 && t_val > 0 && isreal(t_val) && ~isnan(t_val) && ~isinf(t_val)
% Check if the solution is already found (to avoid duplicates)
if isempty(current_fixed_points) || all(abs(current_fixed_points - t_val) > 1e-6)
% Calculate the partial derivative (Jacobian) at the fixed point to check stability
df_dt = (f(t_val + 1e-6, m_val) - f(t_val, m_val)) / 1e-6;
% The eigenvalue in this 1D case is just df_dt
eigenvalue = df_dt;
% Check stability based on the real part of the eigenvalue
is_stable = real(eigenvalue) < 0;
% Store the fixed point and its stability
current_fixed_points = [current_fixed_points, t_val];
current_stability_info = [current_stability_info, is_stable];
end
end
catch
% Skip this initial guess if it causes an error
continue;
end
end
% Store results for current m
fixed_points{i} = current_fixed_points;
stability_info{i} = current_stability_info;
end
% Display results
for i = 1:length(m_values)
m_val = m_values(i);
t_vals = fixed_points{i};
stability_vals = stability_info{i};
for j = 1:length(t_vals)
fprintf('For m = %.2f, t = %.6f: stable = %d\n', m_val, t_vals(j), stability_vals(j));
end
end
% Plot fixed points as a function of m with stability info
figure;
hold on;
for i = 1:length(m_values)
t_vals = fixed_points{i};
stability_vals = stability_info{i};
for j = 1:length(t_vals)
if stability_vals(j)
plot(m_values(i), t_vals(j), 'go'); % Stable points in green
else
plot(m_values(i), t_vals(j), 'ro'); % Unstable points in red
end
end
end
xlabel('m');
ylabel('Fixed Point t');
title('Fixed Points and Stability as a Function of m');
legend('Stable', 'Unstable');
hold off;
21 Commenti
David Goodmanson
il 6 Giu 2024
Hi Sabrina,
I ran the code, and after -- in the second part -- I changed a bunch of capital letters to small letters and funny single quote marks to good single quote marks (I don't know how that part could have run), everything agreed and came out Stable.
Sabrina Garland
il 6 Giu 2024
Sabrina Garland
il 6 Giu 2024
J is always empty in your first code (see above) - thus I guess that this is the reason that the eigenvalue check is always positive.
The command "diff" to differentiate a function is not applicable in numerical computations - at least not for the purpose of your code.
Sabrina Garland
il 6 Giu 2024
Why not using the same method as in your second code ?
J(1,1) = (eq1([k_solution+1e-6, t_solution],m)-eq1([k_solution, t_solution],m))/1e-6;
J(1,2) = (eq1([k_solution, t_solution+1e-6],m)-eq1([k_solution, t_solution],m))/1e-6;
J(2,1) = (eq2([k_solution+1e-6, t_solution],m)-eq2([k_solution, t_solution],m))/1e-6;
J(2,2) = (eq2([k_solution, t_solution+1e-6],m)-eq2([k_solution, t_solution],m))/1e-6;
Sabrina Garland
il 6 Giu 2024
Modificato: Torsten
il 6 Giu 2024
I did and I got entire system as unstable. So, I switched. Here's the code I worked with -
% Define the equations as anonymous functions
eq1 = @(vars, m) (sqrt(vars(1))/sqrt(1-vars(1))) + ((vars(2)*(2*(m^3) + vars(2)^3 - 3*vars(2)*((m+1)^2))^2)/(6*vars(2) - 3*((m-vars(2))^2)) + ((m-vars(2))^2)*(m+2*vars(2))*(m/(3*vars(2)) + 2/3)) / ((m^3 - (vars(2)^2)*(3*m+3) + 2*(vars(2)^3))*(m/(3*vars(2)) - (2*m+vars(2))/(3*(2*vars(2)-((m-vars(2))^2))) + 2/3) - vars(2)*((m-vars(2))^2)*(2*m+vars(2))*((2*m/3) + vars(2)/3));
eq2 = @(vars, m) vars(2) - ((sqrt(vars(1))/sqrt(1-vars(1)))*(((m/vars(2))+2)/3) - (1/3)*(2 + (m/vars(2)) + ((2*m+vars(2))/(((m-vars(2))^2)-2*vars(2))))) / ((sqrt(vars(1))/sqrt(1-vars(1)))*(2*(m^3) - 3*((1+m)^2)*vars(2) + vars(2)^3)/(3*(((m-vars(2))^2)-2*vars(2))) - ((2*m+vars(2))/3));
% Define the range of m values
m_values = linspace(0, 1, 100); % 100 values of m between 0 and 1
% Preallocate arrays to store solutions
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));
stabilities = cell(size(m_values));
% Loop over m values and solve the system of equations for each m
for i = 1:length(m_values)
m = m_values(i);
% System of equations for current m
system_of_equations = @(vars) [eq1(vars, m); eq2(vars, m)];
% Initial guess for [k, t]
initial_guess = [0.75, 1.5];
% Solve the system of equations using fsolve
options = optimoptions('fsolve', 'Display', 'off');
solution = fsolve(system_of_equations, initial_guess, options);
% Store solutions if they meet the criteria
k_solution = solution(1);
t_solution = solution(2);
if k_solution > 0.5 && k_solution < 1 && t_solution > 0
k_solutions(i) = k_solution;
t_solutions(i) = t_solution;
% Compute the Jacobian matrix using finite differences
delta = 1e-6;
J = zeros(2); % Initialize Jacobian matrix
% Calculate the partial derivatives for the Jacobian matrix
for j = 1:2
perturbed_solution = solution;
perturbed_solution(j) = perturbed_solution(j) + delta;
f_perturbed = system_of_equations(perturbed_solution);
f_original = system_of_equations(solution);
J(:, j) = (f_perturbed - f_original) / delta;
end
% Calculate eigenvalues of the Jacobian matrix
eigenvalues = eig(J);
% Analyze stability
if all(real(eigenvalues) < 0)
stabilities{i} = 'Stable';
elseif all(real(eigenvalues) == 0)
% Use center manifold method for stability analysis
% Modify this part with your center manifold method implementation
stabilities{i} = 'Center manifold method';
elseif any(real(eigenvalues) > 0)
stabilities{i} = 'Unstable';
end
else
k_solutions(i) = NaN;
t_solutions(i) = NaN;
stabilities{i} = 'NaN';
end
end
% Display solutions and stabilities
disp('Solutions for each m:');
for i = 1:length(m_values)
fprintf('m = %.2f, k = %.4f, t = %.4f, Stability = %s\n', m_values(i), k_solutions(i), t_solutions(i), stabilities{i});
end
% Plot solutions
figure;
plot(m_values, k_solutions, '-o', 'DisplayName', 'k');
hold on;
plot(m_values, t_solutions, '-o', 'DisplayName', 't');
xlabel('m');
ylabel('Values');
title('Solutions as a Function of m');
legend('k', 't');
hold off;
% Plot stability
figure;
plot(m_values, strcmp(stabilities, 'Stable'), '-o', 'DisplayName', 'Stable');
hold on;
plot(m_values, strcmp(stabilities, 'Unstable'), '-o', 'DisplayName', 'Unstable');
xlabel('m');
ylabel('Stability');
title('Stability as a Function of m');
legend('Stable', 'Unstable');
hold off;
Sabrina Garland
il 6 Giu 2024
A little improvement (not so much relevant in the present case because you only work with two equations):
You should compute
f_original = system_of_equations(solution);
once outside the j-loop and not compute it repeatedly within the j-loop.
And you should check the exitflag from "fsolve" to see if it really converged.
I suggest the following code where I changed the initial guess to the solution for the preceeding m-value:
% Define the equations as anonymous functions
eq1 = @(vars, m) (sqrt(vars(1))/sqrt(1-vars(1))) + ((vars(2)*(2*(m^3) + vars(2)^3 - 3*vars(2)*((m+1)^2))^2)/(6*vars(2) - 3*((m-vars(2))^2)) + ((m-vars(2))^2)*(m+2*vars(2))*(m/(3*vars(2)) + 2/3)) / ((m^3 - (vars(2)^2)*(3*m+3) + 2*(vars(2)^3))*(m/(3*vars(2)) - (2*m+vars(2))/(3*(2*vars(2)-((m-vars(2))^2))) + 2/3) - vars(2)*((m-vars(2))^2)*(2*m+vars(2))*((2*m/3) + vars(2)/3));
eq2 = @(vars, m) vars(2) - ((sqrt(vars(1))/sqrt(1-vars(1)))*(((m/vars(2))+2)/3) - (1/3)*(2 + (m/vars(2)) + ((2*m+vars(2))/(((m-vars(2))^2)-2*vars(2))))) / ((sqrt(vars(1))/sqrt(1-vars(1)))*(2*(m^3) - 3*((1+m)^2)*vars(2) + vars(2)^3)/(3*(((m-vars(2))^2)-2*vars(2))) - ((2*m+vars(2))/3));
% Define the range of m values
m_values = linspace(0, 1, 100); % 100 values of m between 0 and 1
% Preallocate arrays to store solutions
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));
stabilities = cell(size(m_values));
% Initial guess for [k, t]
initial_guess = [0.75, 1.5];
% Loop over m values and solve the system of equations for each m
for i = 1:length(m_values)
m = m_values(i);
% System of equations for current m
system_of_equations = @(vars) [eq1(vars, m); eq2(vars, m)];
% Solve the system of equations using fsolve
options = optimoptions('fsolve', 'Display', 'off');
[solution,~,exitflag] = fsolve(system_of_equations, initial_guess, options);
%exitflag
% Store solutions if they meet the criteria
k_solution = solution(1);
t_solution = solution(2);
if k_solution > 0.5 && k_solution < 1 && t_solution > 0
k_solutions(i) = k_solution;
t_solutions(i) = t_solution;
% Compute the Jacobian matrix using finite differences
delta = 1e-6;
J = zeros(2); % Initialize Jacobian matrix
f_original = system_of_equations(solution);
% Calculate the partial derivatives for the Jacobian matrix
for j = 1:2
perturbed_solution = solution;
perturbed_solution(j) = perturbed_solution(j) + delta;
f_perturbed = system_of_equations(perturbed_solution);
J(:, j) = (f_perturbed - f_original) / delta;
end
% Calculate eigenvalues of the Jacobian matrix
eigenvalues = eig(J);
% Analyze stability
if all(real(eigenvalues) < 0)
stabilities{i} = 'Stable';
elseif all(real(eigenvalues) == 0)
% Use center manifold method for stability analysis
% Modify this part with your center manifold method implementation
stabilities{i} = 'Center manifold method';
elseif any(real(eigenvalues) > 0)
stabilities{i} = 'Unstable';
end
else
k_solutions(i) = NaN;
t_solutions(i) = NaN;
stabilities{i} = 'NaN';
end
% Initial guess for [k, t]
initial_guess = solution;
end
% Display solutions and stabilities
disp('Solutions for each m:');
for i = 1:length(m_values)
fprintf('m = %.2f, k = %.4f, t = %.4f, Stability = %s\n', m_values(i), k_solutions(i), t_solutions(i), stabilities{i});
end
% Plot solutions
figure;
plot(m_values, k_solutions, '-o', 'DisplayName', 'k');
hold on;
plot(m_values, t_solutions, '-o', 'DisplayName', 't');
xlabel('m');
ylabel('Values');
title('Solutions as a Function of m');
legend('k', 't');
hold off;
% Plot stability
figure;
plot(m_values, strcmp(stabilities, 'Stable'), '-o', 'DisplayName', 'Stable');
hold on;
plot(m_values, strcmp(stabilities, 'Unstable'), '-o', 'DisplayName', 'Unstable');
xlabel('m');
ylabel('Stability');
title('Stability as a Function of m');
legend('Stable', 'Unstable');
hold off;
Sabrina Garland
il 6 Giu 2024
Modificato: Torsten
il 6 Giu 2024
Thanks for suggestion. I am still stuck.
% Define the equations as anonymous functions
eq1 = @(vars, m) (sqrt(vars(1))/sqrt(1-vars(1))) + ((vars(2)*(2*(m^3) + vars(2)^3 - 3*vars(2)*((m+1)^2))^2)/(6*vars(2) - 3*((m-vars(2))^2)) + ((m-vars(2))^2)*(m+2*vars(2))*(m/(3*vars(2)) + 2/3)) / ((m^3 - (vars(2)^2)*(3*m+3) + 2*(vars(2)^3))*(m/(3*vars(2)) - (2*m+vars(2))/(3*(2*vars(2)-((m-vars(2))^2))) + 2/3) - vars(2)*((m-vars(2))^2)*(2*m+vars(2))*((2*m/3) + vars(2)/3));
eq2 = @(vars, m) vars(2) - ((sqrt(vars(1))/sqrt(1-vars(1)))*(((m/vars(2))+2)/3) - (1/3)*(2 + (m/vars(2)) + ((2*m+vars(2))/(((m-vars(2))^2)-2*vars(2))))) / ((sqrt(vars(1))/sqrt(1-vars(1)))*(2*(m^3) - 3*((1+m)^2)*vars(2) + vars(2)^3)/(3*(((m-vars(2))^2)-2*vars(2))) - ((2*m+vars(2))/3));
% Define the range of m values
m_values = linspace(0, 1, 100); % 100 values of m between 0 and 1
% Preallocate arrays to store solutions
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));
stabilities = cell(size(m_values));
% Loop over m values and solve the system of equations for each m
for i = 1:length(m_values)
m = m_values(i);
% System of equations for current m
system_of_equations = @(vars) [eq1(vars, m); eq2(vars, m)];
% Initial guess for [k, t]
initial_guess = [0.75, 1.5];
% Solve the system of equations using fsolve
options = optimoptions('fsolve', 'Display', 'off');
[solution, ~, exitflag] = fsolve(system_of_equations, initial_guess, options);
% Check if fsolve converged
if exitflag > 0
% Store solutions if they meet the criteria
k_solution = solution(1);
t_solution = solution(2);
if k_solution > 0.5 && k_solution < 1 && t_solution > 0
k_solutions(i) = k_solution;
t_solutions(i) = t_solution;
% Compute the Jacobian matrix using finite differences
delta = 1e-6;
J = zeros(2); % Initialize Jacobian matrix
% Compute f_original once outside the loop
f_original = system_of_equations(solution);
% Calculate the partial derivatives for the Jacobian matrix
for j = 1:2
perturbed_solution = solution;
perturbed_solution(j) = perturbed_solution(j) + delta;
f_perturbed = system_of_equations(perturbed_solution);
J(:, j) = (f_perturbed - f_original) / delta;
end
% Calculate eigenvalues of the Jacobian matrix
eigenvalues = eig(J);
% Analyze stability
if all(real(eigenvalues) < 0)
stabilities{i} = 'Stable';
elseif all(real(eigenvalues) == 0)
% Use center manifold method for stability analysis
% Modify this part with your center manifold method implementation
stabilities{i} = 'Center manifold method';
elseif any(real(eigenvalues) > 0)
stabilities{i} = 'Unstable';
end
else
k_solutions(i) = NaN;
t_solutions(i) = NaN;
stabilities{i} = 'NaN';
end
else
k_solutions(i) = NaN;
t_solutions(i) = NaN;
stabilities{i} = 'NaN';
end
end
% Display solutions and stabilities
disp('Solutions for each m:');
for i = 1:length(m_values)
fprintf('m = %.2f, k = %.4f, t = %.4f, Stability = %s\n', m_values(i), k_solutions(i), t_solutions(i), stabilities{i});
end
% Plot solutions
figure;
plot(m_values, k_solutions, '-o', 'DisplayName', 'k');
hold on;
plot(m_values, t_solutions, '-o', 'DisplayName', 't');
xlabel('m');
ylabel('Values');
title('Solutions as a Function of m');
legend('k', 't');
hold off;
% Plot stability
figure;
plot(m_values, strcmp(stabilities, 'Stable'), '-o', 'DisplayName', 'Stable');
hold on;
plot(m_values, strcmp(stabilities, 'Unstable'), '-o', 'DisplayName', 'Unstable');
xlabel('m');
ylabel('Stability');
title('Stability as a Function of m');
legend('Stable', 'Unstable');
hold off;
Torsten
il 6 Giu 2024
And where are you stuck ?
Sabrina Garland
il 6 Giu 2024
Sabrina Garland
il 6 Giu 2024
Modificato: Sabrina Garland
il 6 Giu 2024
Given the nonlinear system
% Define the equations as anonymous functions
eq1 = @(vars, m) (sqrt(vars(1))/sqrt(1-vars(1))) + ((vars(2)*(2*(m^3) + vars(2)^3 - 3*vars(2)*((m+1)^2))^2)/(6*vars(2) - 3*((m-vars(2))^2)) + ((m-vars(2))^2)*(m+2*vars(2))*(m/(3*vars(2)) + 2/3)) / ((m^3 - (vars(2)^2)*(3*m+3) + 2*(vars(2)^3))*(m/(3*vars(2)) - (2*m+vars(2))/(3*(2*vars(2)-((m-vars(2))^2))) + 2/3) - vars(2)*((m-vars(2))^2)*(2*m+vars(2))*((2*m/3) + vars(2)/3));
eq2 = @(vars, m) vars(2) - ((sqrt(vars(1))/sqrt(1-vars(1)))*(((m/vars(2))+2)/3) - (1/3)*(2 + (m/vars(2)) + ((2*m+vars(2))/(((m-vars(2))^2)-2*vars(2))))) / ((sqrt(vars(1))/sqrt(1-vars(1)))*(2*(m^3) - 3*((1+m)^2)*vars(2) + vars(2)^3)/(3*(((m-vars(2))^2)-2*vars(2))) - ((2*m+vars(2))/3));
you only get unstable solutions for 0 <= m <= 1. That's understandable from your code.
Maybe you could explain what you do in your second approach and how this is related to the first investigation for stability. Since your function f has changed, I can't understand why you think both approaches will give the same result, especially since you keep k constant as k = 0.6945.
Steven Lord
il 6 Giu 2024
In the cases where you expected the code to report stability but it reported an unstable result, what are the values of the positive eigenvalues? Are they much greater than zero or are they essentially zero but because of numerical noise are just barely greater than zero?
Sabrina Garland
il 6 Giu 2024
Sabrina Garland
il 6 Giu 2024
Sabrina Garland
il 6 Giu 2024
Torsten
il 7 Giu 2024
As I could see, eigenvalues for t are near zero but eigenvalues for k are quite large.
The equations are not separable - there are no eigenvalues for t and eigenvalues for k. So you cannot say that the system is stable with respect to t and unstable with respect to k or something similar.
Sabrina Garland
il 7 Giu 2024
Risposte (0)
Categorie
Scopri di più su Choose a Solver in Centro assistenza e File Exchange
il 6 Giu 2024
il 7 Giu 2024
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
