Mismatched Stability check.
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 = %sn’, 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 = %dn’, 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;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 = %sn’, 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 = %dn’, 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; 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 = %sn’, 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 = %dn’, 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; matlab, optimization MATLAB Answers — New Questions
