Vous suivez désormais cette question
- Les mises à jour seront visibles dans votre flux de contenu suivi.
- Selon vos préférences en matière de communication il est possible que vous receviez des e-mails.
Mismatched Stability check.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710301/image.png)
21 commentaires
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710351/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710356/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710601/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710606/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710371/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1710376/image.png)
Here's the code checking the eigenvalues. As I could see, eigenvalues for t are near zero but eigenvalues for k are quite large.
Should I take different perturbation for k and t?
% 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);
% Debug: print eigenvalues for inspection fprintf('m = %.2f, k = %.4f, t = %.4f, Eigenvalues: %.6f, %.6f\n', m, k_solution, t_solution, eigenvalues(1), eigenvalues(2));
% 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;
Réponses (0)
Voir également
Catégories
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!Une erreur s'est produite
Impossible de terminer l’action en raison de modifications de la page. Rechargez la page pour voir sa mise à jour.
Sélectionner un site web
Choisissez un site web pour accéder au contenu traduit dans votre langue (lorsqu'il est disponible) et voir les événements et les offres locales. D’après votre position, nous vous recommandons de sélectionner la région suivante : .
Vous pouvez également sélectionner un site web dans la liste suivante :
Comment optimiser les performances du site
Pour optimiser les performances du site, sélectionnez la région Chine (en chinois ou en anglais). Les sites de MathWorks pour les autres pays ne sont pas optimisés pour les visites provenant de votre région.
Amériques
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asie-Pacifique
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)