clc;
close all;
%% parameter
R0=0.8; a0=2923; a1=-628; a2=402.1; Cth=1.324e-13; Tth=5.953e-6; Tamb=298;
CC=100;
Vm=0.30032043; T=701.5886849;
%求雅可比矩阵
%% Jacobian
for i=1:length(CC)
C=CC(i,1);
A11=(-1./C.*10e9)*(1-Vm.*(a1./(2.*T.*Vm.^0.5)+a2./T))./(R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T));
A12=(-1./C.*10e9).*(Vm.*(a0+a1.*Vm.^0.5+a2.*Vm))./(T.^2.*(R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T)));
B11=(2.*Vm-Vm.^2.*(a1./(2.*T.*Vm.^0.5)+a2./T))./(Cth.*(R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T)));
B12=((Vm.^2.*(a0+a1.*Vm.^0.5+a2.*Vm))./(T.^2.*R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T))-Tth)./Cth;
A=[A11 A12;
B11 B12];
[V,D]=eig(A);%求矩阵的特征值特征向量
X=diag(D);
xx(i)=X(1,:);
yy(i)=X(2,:);
XX=xx';
YY=yy';
Re1=real(XX);
Im1=imag(XX);
Re2=real(YY);
Im2=imag(YY);
end
根据你提供的代码,我看到 xx 和 yy 是在循环内部定义的数组,而 XX 和 YY 是在循环结束后定义的数组,且是通过将 xx 和 yy 转置得到的。因此,在你的代码中,只能得到最后一次循环的结果。如果你想要得到所有循环的结果,需要将 XX 和 YY 改为矩阵形式的变量,并在循环中不断将每次的结果添加到这些矩阵中。
以下是我帮你修改后的代码示例,可以试试:
clc;
close all;
%% parameter
R0=0.8; a0=2923; a1=-628; a2=402.1; Cth=1.324e-13; Tth=5.953e-6; Tamb=298;
CC=100;
Vm=0.30032043; T=701.5886849;
%% Jacobian
num_loop = length(CC);
A_matrix = zeros(num_loop, 2, 2); % 定义矩阵形式的变量来存储 A 矩阵
eig_values = zeros(num_loop, 2);
for i = 1:num_loop
C = CC(i,1);
A11 = (-1./C.*10e9)*(1-Vm.*(a1./(2.*T.*Vm.^0.5)+a2./T))./(R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T));
A12 = (-1./C.*10e9).*(Vm.*(a0+a1.*Vm.^0.5+a2.*Vm))./(T.^2.*(R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T)));
B11 = (2.*Vm-Vm.^2.*(a1./(2.*T.*Vm.^0.5)+a2./T))./(Cth.*(R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T)));
B12 = ((Vm.^2.*(a0+a1.*Vm.^0.5+a2.*Vm))./(T.^2.*R0.*exp((a0+a1.*Vm.^0.5+a2.*Vm)./T))-Tth)./Cth;
A = [A11 A12; B11 B12];
A_matrix(i,:,:) = A; % 将每次计算得到的 A 存储到矩阵中
[V, D] = eig(A);
eig_values(i,:) = diag(D);
end
Re1 = real(eig_values(:,1));
Im1 = imag(eig_values(:,1));
Re2 = real(eig_values(:,2));
Im2 = imag(eig_values(:,2));
% 绘制结果的代码
figure;
subplot(2,2,1); scatter(Re1, Im1, 10, 'filled'); xlabel('Real'); ylabel('Imaginary'); title('Eigenvalue 1');
subplot(2,2,2); scatter(Re2, Im2, 10, 'filled'); xlabel('Real'); ylabel('Imaginary'); title('Eigenvalue 2');
subplot(2,2,3); scatter(Re1, Re2, 10, 'filled'); xlabel('Eigenvalue
感谢回答!还有一点问题,关于for循环这个问题已成功解决!
这个程序我的想法是CC变化求解雅可比矩阵特征值,λ1、λ2,然后再得到它们对应的虚部、实部数据,最后绘制出横坐标为特征值实部,纵坐标为特征值虚部得,随着CC得值变化得曲线,类似下图。
但是特征值得应用部分好像不正确,可以麻烦你帮我改改吗
见正文回复
追问你好,非常感谢前面的回答!再问一下,我更换了方程组和参数,现在的问题是[V,D]=eig(A)这个矩阵不是方阵,但是不知道怎么改,A很显然就是方阵呀,