回复 2 # yuexiaoming 的帖子
给你一个自编的判断解是什么类型的奇点程序(程序的内部算法源于刘延柱老师的非线性振动那本书的第二章)
%求解系统的极点和并判断其类型
%求解的系统类型:
% dx1=f(x1,x2);
% dx2=g(x1,x2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function CC=critical_points()
clc
clear all
syms x1 x2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[F,V Cri]=aa(); %奇点
%F是两个函数组成的矩阵,V是函数中得变量V=[x1 x2];Cri为奇点
J=jacobin(F,V);
for k=1:1:size(Cri,1)
A=subs(J,{x1,x2},{Cri(k,1),Cri(k,2)});
A=double(A);
p=trace(A);
q=det(A) ; %节点类型取决于A的特征值结构
w=p^2-4*q;
%奇点类型判断%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if w>=0 %w>0,r1和r2 是不等实根,w=0是重根
if q>0 %R1和R2是同号
if p<=0
Cr=[Cri(k,1),Cri(k,2)]
disp('稳定节点')
elseif p>0
Cr=[Cri(k,1),Cri(k,2)]
disp('不稳定节点')
end
elseif q<0
Cr=[Cri(k,1),Cri(k,2)]
disp('鞍点')
end
elseif w<0 %共轭复根
if p==0
Cr=[Cri(k,1),Cri(k,2)]
disp('中点')
else p~=0
if p<=0
Cr=[Cri(k,1),Cri(k,2)]
disp('稳定焦点')
elseif p>0
Cr=[Cri(k,1),Cri(k,2)]
disp('不稳定焦点')
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
%%%%%%%%%%%不同的函数需在此设置
function [F,V Cri]=aa()
%F是两个函数组成的矩阵,V是函数中得变量V=[x1 x2];Cri为奇点
syms x1 x2
F=[x1*(1-x1-x2);0.25*x2*(2-3*x1-x2)];%两个函数就是两行
V=[x1,x2];
[x1,x2]=solve('x1*(1-x1-x2)=0','0.25*x2*(2-3*x1-x2)= 0'); %函数,不同的需更改
Cri=[x1,x2];
end
|