% %% ========== 微分方程相图绘制大全 ==========
% % 包含一维相线、二维线性系统、多种经典非线性系统
% % 直接运行本脚本即可依次显示所有图形
% 
% clear; clc; close all;
% 
% %% ---------- 1. 一维相线：Logistic 方程 dx/dt = x(1-x) ----------
% f1 = @(x) x.*(1-x);
% x_vals = linspace(-0.5, 1.5, 200);
% dx_vals = f1(x_vals);
% 
% figure('Name', '1D 相线：Logistic 方程', 'NumberTitle', 'off', 'Position', [50 550 700 150]);
% plot(x_vals, zeros(size(x_vals)), 'k-', 'LineWidth', 1.5); hold on;
% % 平衡点
% x_eq = [0, 1];
% for i = 1:length(x_eq)
%     % 判断稳定性：导数符号变化
%     if f1(x_eq(i)-0.01) > 0 && f1(x_eq(i)+0.01) < 0
%         plot(x_eq(i), 0, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 8); % 稳定
%     elseif f1(x_eq(i)-0.01) < 0 && f1(x_eq(i)+0.01) > 0
%         plot(x_eq(i), 0, 'ko', 'MarkerFaceColor', 'w', 'MarkerSize', 8); % 不稳定
%     else
%         plot(x_eq(i), 0, 'ko', 'MarkerFaceColor', [0.5 0.5 0.5], 'MarkerSize', 8); % 半稳定
%     end
% end
% % 箭头
% for xi = linspace(-0.3, 1.3, 15)
%     dir = sign(f1(xi));
%     if dir > 0
%         quiver(xi, 0, 0.04, 0, 'b', 'MaxHeadSize', 0.6, 'LineWidth', 1);
%     elseif dir < 0
%         quiver(xi, 0, -0.04, 0, 'b', 'MaxHeadSize', 0.6, 'LineWidth', 1);
%     end
% end
% xlabel('x'); ylabel(''); title('dx/dt = x(1-x)   ●稳定   ○不稳定');
% ylim([-0.2 0.2]); grid on; set(gca, 'YTick', []);
% 
% 
% %% ---------- 2. 二维线性系统：6 种典型相图 ----------
% % 定义系统矩阵和名称
% linear_systems = {
%     struct('A', [-1 0; 0 -2], 'name', '稳定结点');
%     struct('A', [1 0; 0 2], 'name', '不稳定结点');
%     struct('A', [1 1; 0 -2], 'name', '鞍点');
%     struct('A', [-0.5 1; -1 -0.5], 'name', '稳定焦点');
%     struct('A', [0.5 1; -1 0.5], 'name', '不稳定焦点');
%     struct('A', [0 1; -1 0], 'name', '中心');
% };
% 
% figure('Name', '二维线性系统相图集合', 'NumberTitle', 'off', 'Position', [100 50 1200 700]);
% for k = 1:6
%     subplot(2, 3, k);
%     A = linear_systems{k}.A;
%     % 网格
%     [X, Y] = meshgrid(-2:0.25:2, -2:0.25:2);
%     U = A(1,1)*X + A(1,2)*Y;
%     V = A(2,1)*X + A(2,2)*Y;
%     % 归一化箭头长度以便显示
%     L = sqrt(U.^2 + V.^2);
%     U = U ./ L; V = V ./ L;
%     quiver(X, Y, U, V, 0.4, 'b'); hold on;
% 
%     % 画轨线
%     tspan = [0 8];
%     if k == 6  % 中心需要反向时间也画
%         tspan_back = [0 -8];
%     end
%     init_grid = -2:0.8:2;
%     for x0 = init_grid
%         for y0 = init_grid
%             if x0==0 && y0==0, continue; end
%             [~, traj] = ode45(@(t,z) A*z, tspan, [x0; y0]);
%             plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1.2);
%             if k == 6
%                 [~, traj] = ode45(@(t,z) A*z, tspan_back, [x0; y0]);
%                 plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1.2);
%             end
%         end
%     end
%     % 画特征向量方向（如果特征值为实数）
%     [Vect, D] = eig(A);
%     if isreal(D)
%         for j = 1:2
%             v = Vect(:,j);
%             line([-3*v(1) 3*v(1)], [-3*v(2) 3*v(2)], 'Color', 'g', 'LineWidth', 1.5, 'LineStyle', '--');
%         end
%     end
%     % 标记平衡点
%     plot(0, 0, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 6);
%     xlabel('x'); ylabel('y'); title(linear_systems{k}.name);
%     axis equal; xlim([-2.5 2.5]); ylim([-2.5 2.5]); grid on;
% end
% 
% 
% %% ---------- 3. 非线性系统：Lotka-Volterra 捕食者-猎物模型 ----------
% % dx/dt = x*(1 - y)    (猎物)
% % dy/dt = y*(x - 1)    (捕食者)
% lv_sys = @(t, z) [z(1)*(1 - z(2)); z(2)*(z(1) - 1)];
% 
% figure('Name', 'Lotka-Volterra 捕食者-猎物模型', 'NumberTitle', 'off', 'Position', [200 50 500 450]);
% [X, Y] = meshgrid(0:0.2:2.5, 0:0.2:2.5);
% U = X.*(1 - Y);
% V = Y.*(X - 1);
% L = sqrt(U.^2 + V.^2) + 1e-8;
% quiver(X, Y, U./L, V./L, 0.5, 'b'); hold on;
% 
% % 零倾线
% fimplicit(@(x,y) x.*(1-y), [0 2.5 0 2.5], 'g--', 'LineWidth', 1.5); % x-nullcline
% fimplicit(@(x,y) y.*(x-1), [0 2.5 0 2.5], 'm--', 'LineWidth', 1.5); % y-nullcline
% 
% % 平衡点
% plot(1, 1, 'ko', 'MarkerFaceColor', 'k', 'MarkerSize', 8); % 中心（非线性下是中心）
% plot(0, 0, 'ko', 'MarkerFaceColor', 'w', 'MarkerSize', 8); % 鞍点
% 
% % 轨线
% for x0 = 0.2:0.4:2.2
%     for y0 = 0.2:0.4:2.2
%         if x0==0 && y0==0, continue; end
%         [~, traj] = ode45(lv_sys, [0 20], [x0; y0]);
%         plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1);
%     end
% end
% xlabel('猎物 x'); ylabel('捕食者 y'); title('Lotka-Volterra 模型 (中心)');
% legend({'向量场','x-零倾线','y-零倾线','轨线'}, 'Location', 'best');
% axis equal; xlim([0 2.5]); ylim([0 2.5]); grid on;
% 
% 
% %% ---------- 4. 非线性系统：竞争模型 ----------
% % dx/dt = x*(2 - x - y)
% % dy/dt = y*(2 - 2*x - y)
% comp_sys = @(t, z) [z(1)*(2 - z(1) - z(2)); z(2)*(2 - 2*z(1) - z(2))];
% 
% figure('Name', '竞争模型 (两个物种)', 'NumberTitle', 'off', 'Position', [720 50 500 450]);
% [X, Y] = meshgrid(0:0.2:2.5, 0:0.2:2.5);
% U = X.*(2 - X - Y);
% V = Y.*(2 - 2*X - Y);
% L = sqrt(U.^2 + V.^2) + 1e-8;
% quiver(X, Y, U./L, V./L, 0.5, 'b'); hold on;
% 
% % 零倾线
% fimplicit(@(x,y) x.*(2-x-y), [0 2.5 0 2.5], 'g--', 'LineWidth', 1.5);
% fimplicit(@(x,y) y.*(2-2*x-y), [0 2.5 0 2.5], 'm--', 'LineWidth', 1.5);
% 
% % 平衡点
% plot(0,0,'ko','MarkerFaceColor','w','MarkerSize',8);     % 不稳定结点
% plot(2,0,'ko','MarkerFaceColor','r','MarkerSize',8);     % 稳定结点（物种1胜）
% plot(0,2,'ko','MarkerFaceColor','r','MarkerSize',8);     % 稳定结点（物种2胜）
% plot(2/3, 4/3,'ko','MarkerFaceColor','w','MarkerSize',8);% 鞍点
% 
% % 轨线
% for x0 = 0.2:0.5:2.2
%     for y0 = 0.2:0.5:2.2
%         if x0==0 && y0==0, continue; end
%         [~, traj] = ode45(comp_sys, [0 15], [x0; y0]);
%         plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1);
%     end
% end
% xlabel('物种 1'); ylabel('物种 2'); title('竞争模型：双稳定结点 + 鞍点');
% legend({'向量场','x-零倾线','y-零倾线','轨线'});
% axis equal; xlim([0 2.5]); ylim([0 2.5]); grid on;
% 
% 
% %% ---------- 5. 非线性系统：Duffing 振子 (双势阱) ----------
% % 等价于一阶系统：
% % dx/dt = y
% % dy/dt = x - x^3 - delta*y   (取 delta = 0.2)
% delta = 0.2;
% duffing = @(t, z) [z(2); z(1) - z(1)^3 - delta*z(2)];
% 
% figure('Name', 'Duffing 振子 (双势阱)', 'NumberTitle', 'off', 'Position', [300 50 500 450]);
% [X, Y] = meshgrid(-2:0.2:2, -1.5:0.2:1.5);
% U = Y;
% V = X - X.^3 - delta*Y;
% L = sqrt(U.^2 + V.^2) + 1e-8;
% quiver(X, Y, U./L, V./L, 0.5, 'b'); hold on;
% 
% % 零倾线
% fimplicit(@(x,y) y, [-2 2 -1.5 1.5], 'g--', 'LineWidth', 1.5);          % x-nullcline
% fimplicit(@(x,y) x - x.^3 - delta*y, [-2 2 -1.5 1.5], 'm--', 'LineWidth', 1.5); % y-nullcline
% 
% % 平衡点
% plot(-1,0,'ro','MarkerFaceColor','r','MarkerSize',8); % 稳定焦点（左势阱）
% plot(1,0,'ro','MarkerFaceColor','r','MarkerSize',8);  % 稳定焦点（右势阱）
% plot(0,0,'ko','MarkerFaceColor','w','MarkerSize',8);  % 鞍点（势垒顶部）
% 
% % 轨线
% for x0 = -1.8:0.6:1.8
%     for y0 = -1.2:0.6:1.2
%         [~, traj] = ode45(duffing, [0 20], [x0; y0]);
%         plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1);
%         [~, traj] = ode45(duffing, [0 -10], [x0; y0]);
%         plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1);
%     end
% end
% xlabel('x'); ylabel('y'); title('Duffing 振子：双稳定焦点 + 鞍点');
% legend({'向量场','x-零倾线','y-零倾线','轨线'});
% axis equal; xlim([-2 2]); ylim([-1.5 1.5]); grid on;
% 
% 
% %% ---------- 6. 非线性系统：Van der Pol 振子 (极限环) ----------
% % dx/dt = y
% % dy/dt = mu*(1 - x^2)*y - x   (取 mu = 1)
% mu = 1;
% vdp = @(t, z) [z(2); mu*(1 - z(1)^2)*z(2) - z(1)];
% 
% figure('Name', 'Van der Pol 振子 (极限环)', 'NumberTitle', 'off', 'Position', [820 50 500 450]);
% [X, Y] = meshgrid(-3:0.25:3, -3:0.25:3);
% U = Y;
% V = mu*(1 - X.^2).*Y - X;
% L = sqrt(U.^2 + V.^2) + 1e-8;
% quiver(X, Y, U./L, V./L, 0.5, 'b'); hold on;
% 
% % 零倾线
% fimplicit(@(x,y) y, [-3 3 -3 3], 'g--', 'LineWidth', 1.5);
% fimplicit(@(x,y) mu*(1-x.^2).*y - x, [-3 3 -3 3], 'm--', 'LineWidth', 1.5);
% 
% % 平衡点（原点是不稳定焦点）
% plot(0,0,'ko','MarkerFaceColor','w','MarkerSize',8);
% 
% % 轨线：从内部和外部出发，最终都趋向极限环
% init_pos = [0.5 0; 2.5 0; -0.5 0; -2.5 0; 0 0.5; 0 2.5];
% for i = 1:size(init_pos,1)
%     [~, traj] = ode45(vdp, [0 20], init_pos(i,:));
%     plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1.5);
% end
% % 为了更清晰，多画几条
% for r = 0.3:0.6:2.7
%     for theta = 0:pi/4:2*pi
%         x0 = r*cos(theta); y0 = r*sin(theta);
%         [~, traj] = ode45(vdp, [0 20], [x0; y0]);
%         plot(traj(:,1), traj(:,2), 'r', 'LineWidth', 1);
%     end
% end
% xlabel('x'); ylabel('y'); title('Van der Pol 振子：不稳定焦点 + 稳定极限环');
% legend({'向量场','x-零倾线','y-零倾线','轨线'});
% axis equal; xlim([-3 3]); ylim([-3 3]); grid on;
% 
% disp('所有相图绘制完成！');
% 
% % 线性系统相图绘制
% % dx/dt = a*x + b*y
% % dy/dt = c*x + d*y
% clear; clc; close all;

%% ========== 用户参数设置 ==========
% 修改这里的 a,b,c,d 以改变系统
a = 1;   b = -1;  % 示例：中心 (特征值 ±i)
c = 0;   d = 1;
% 其他典型示例（取消注释使用）：
% a = -1; b = 0;  c = 0; d = -1;   % 稳定结点
% a = 1;  b = 1;  c = 4; d = 1;    % 鞍点
% a = -1; b = -2; c = 1; d = -1;   % 稳定焦点

% 绘图范围
x_range = [-2, 2];
y_range = [-2, 2];

% 向量场网格密度
N = 20;

% 积分时间 (可根据系统调整，发散过快请减小 t_end)
t_start = 0;
t_end   = 5;

% 初始点集合 (可自行增删)
init_points = [ 1,  1;
               -1, -1;
                1, -1;
               -1,  1;
                1.5, 0;
                0,  1.5;
               -1.5, 0;
                0, -1.5 ];
%% ==================================

% 创建网格
x = linspace(x_range(1), x_range(2), N);
y = linspace(y_range(1), y_range(2), N);
[X, Y] = meshgrid(x, y);

% 计算向量场分量
U = a*X + b*Y;
V = c*X + d*Y;

% 绘制向量场
figure;
quiver(X, Y, U, V, 'r', 'AutoScale', 'on', 'AutoScaleFactor', 0.3);
hold on;

% 绘制零倾线 (dx/dt=0 和 dy/dt=0)
% dx/dt = 0 : a*x + b*y = 0
if abs(b) > 1e-6
    x_null = x_range;
    y_null1 = -a/b * x_null;
    plot(x_null, y_null1, 'g--', 'LineWidth', 1);
elseif abs(a) > 1e-6
    plot([0 0], y_range, 'g--', 'LineWidth', 1);
end
% dy/dt = 0 : c*x + d*y = 0
if abs(d) > 1e-6
    x_null = x_range;
    y_null2 = -c/d * x_null;
    plot(x_null, y_null2, 'b--', 'LineWidth', 1);
elseif abs(c) > 1e-6
    plot(x_range, [0 0], 'b--', 'LineWidth', 1);
end

% 定义微分方程
f = @(t, y) [a*y(1) + b*y(2); c*y(1) + d*y(2)];

% 积分求解相轨迹
tspan = [t_start, t_end];
for i = 1:size(init_points, 1)
    y0 = init_points(i, :)';
    [~, traj] = ode45(f, tspan, y0);
    plot(traj(:,1), traj(:,2), 'LineWidth', 1.5);  % 相轨迹
    plot(y0(1), y0(2), 'ko', 'MarkerFaceColor', 'k'); % 初始点
end

% 标记平衡点 (原点)
plot(0, 0, 'ro', 'MarkerSize', 8, 'MarkerFaceColor', 'r');

% 图形修饰
xlabel('x');
ylabel('y');
title(sprintf('相图: dx/dt = %g x %+g y,  dy/dt = %g x %+g y', a, b, c, d));
grid on;
axis equal;
xlim(x_range);
ylim(y_range);
legend('向量场', 'dx/dt=0', 'dy/dt=0', '相轨迹', '初始点', '平衡点', ...
       'Location', 'best');
hold off;

% 输出特征值信息 (辅助分析)
A = [a, b; c, d];
eigvals = eig(A);
fprintf('系统矩阵的特征值: λ1 = %.4f, λ2 = %.4f\n', eigvals(1), eigvals(2));