reach to phase portrait in a 4*4 matrix and plotting all of them in one figure (if possible), otherwise plotting them 2 by 2 versus to each other.
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi guys. I wanna draw phase portarait of a matrix 4*4 that is solved for each of its 2 states based on each other. when I write this code, it gives me this error:
Error in
@(t,V)[(3.7*10^-15*V(1)+7.6*10^-6*V(2));(-4*10^-16*V(1)-6.5*10^-5*V(2));(2.37*10^-16*V(1)-0.142*V(3));(-8.15*10^-14*V(1)-0.5535*V(2)+18.21*V(3)+0.0833*V(4))]
Error in test111 (line 70)
Yprime = f1(t,[x(j); y(j)]);
And also why when I change the initial conditions to ones which I comment, it does not give me the same solution after 2500 iteration? Actually it is the operation point.(x1(1)=57.5,x2(1)=8.5 x3(1)=0.051, x4(1)=4.8) and I expect that I reach around this states whatever initial condition that I start.
Thanks
clc;
clear all;
close all;
x1(1)=57.5;
x2(1)=8.5;
x3(1)=0.051;
x4(1)=4.8;
% x1(1)=80;
% x2(1)=10;
% x3(1)=1;
% x4(1)=8;
deltat=10^-3;
A1=[3.4*10^-15 7.6*10^-6 0 0;-4.08*10^-16 -6.55*10^-5 0 0;2.37*10^-16 0.00059 -0.14 0;-8.15*10^-14 -0.05 18.216 0.083];
A2=[9.04*10^-7 -1.48*10^-6 0 0;-1.97*10^-7 -1*10^-5 0 0;9.24*10^-8 0.000458 -0.12 0;-8.55*10^-5 -0.104 65.217 -0.083];
B1=[0.0015 -0.0015 -6.9*10^-12;-5.9*10^-5 -9*10^2 5.9*10^-11;3.46*10^-5 5.2*10^-5 3.24*10^-11;-0.0167 0.0239 -1.0733*10^-9];
B2=[0.0014 -0.0013 -1.9*10^-11;-2.9*10^-5 -0.00011 6.1*10^-11;1.36*10^-5 5.33*10^-5 2.15*10^-11;-0.01 0.042 -4.47*10^-8];
u=[32;32;80];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
P=[161469.069263159,8918967.27830892,-576.987918162128,1.43011607072035;8918967.27830892,494939491.276476,-3069582.33611650,-8755.93015894571;-576.987918162128,-3069582.33611650,708027606.278748,206347.846527781;1.43011607072035,-8755.93015894571,206347.846527781,4552.55710715823];
for i=1:2500
% x1(i+1)=(3.7*10^-15*x1(i)+7.6*10^-6*x2(i)+0.554*10^-9)*deltat+x1(i);
% x2(i+1)=(-4*10^-16*x1(i)-6.5*10^-5*x2(i)-2.8)*deltat+x2(i);
% x3(i+1)=(2.37*10^-16*x1(i)-0.142*x3(i)-0.16*10^-2+0.275*10^-2)*deltat+x3(i);
% x4(i+1)=(-8.15*10^-14*x1(i)-0.5535*x2(i)+18.21*x3(i)+0.0833*x4(i)+0.2)*deltat+x4(i);
% V1(i)=0.5*(x1(i).^2+x2(i).^2+x3(i).^2+x4(i).^2);
x1(i+1)=(9.04*10^-7*x1(i)-1.48*10^-6*x2(i)+0.0032)*deltat+x1(i);
x2(i+1)=(-1.97*10^-7*x1(i)-1.008*10^-5*x2(i)-0.0044)*deltat+x2(i);
x3(i+1)=(9.2*10^-8*x1(i)+0.000458*x2(i)-0.127*x3(i)+0.0021)*deltat+x3(i);
x4(i+1)=(-8.55*10^-5*x1(i)-0.1*x2(i)+65.217*x3(i)-0.0833*x4(i)+1.024)*deltat+x4(i);
V2(i)=[x1(i) x2(i) x3(i) x4(i)]*P*[x1(i) x2(i) x3(i) x4(i)]';
%V2(i)=x1(i).^2+300000000*x2(i).^2+200000000*x3(i).^2+0.00000001*x4(i).^2;
end
for i=1:2500
figure(1)
plot(x1(1,:),x2(1,:))
hold on;
plot(x1(1,1),x2(1,1),'bo') % starting point
hold on;
plot(x1(1,end),x2(1,end),'ks') % ending point
hold on;
figure(2);
plot(i,x1(1,i),'*-')
hold on;
% figure(3)
% plot(i,V1(1,i),'*-')
% hold on
figure(4)
plot(i,V2(1,i),'*-')
hold on;
end
hold off;
%---------------------- Vector field & phase portrait
f1=@(t,V)[(3.7*10^-15*V(1)+7.6*10^-6*V(2));(-4*10^-16*V(1)-6.5*10^-5*V(2));(2.37*10^-16*V(1)-0.142*V(3));(-8.15*10^-14*V(1)-0.5535*V(2)+18.21*V(3)+0.0833*V(4))];
f2=@(t,V)[(9.04*10^-7*V(1)-1.48*10^-6*V(2));(-1.97*10^-7*V(1)-1.008*10^-5*V(2));(9.2*10^-8*V(1)+0.000458*V(2)-0.127*V(3));(-8.55*10^-5*V(1)-0.1*V(2)+65.217*V(3)-0.0833*V(4))];
v1 = linspace(-10,10,50);
v2 = linspace(-10,10,50);
v3 = linspace(-10,10,50);
v4 = linspace(-10,10,50);
[x,y] = meshgrid(v1,v2);
size(x)
size(y)
u = zeros(size(x));
v = zeros(size(y));
s = zeros(size(x));
z = zeros(size(y));
t=0; % we want the derivatives at each point at t=0, i.e. the starting time
for j = 1:numel(x)
Yprime = f1(t,[x(j); y(j)]);
Yprim = f2(t,[x(j); y(j)]);
u(j) = Yprime(1);
v(j) = Yprime(2);
s(j)= Yprim(1);
z(j) = Yprim(2);
end
figure(5)
quiver(x,y,u,v,'r');
figure(6)
quiver(x,y,s,z,'b');
xlabel('x1')
ylabel('x2')
legend('vector field')
axis tight equal;
for y10=[0 1]
for y20=[0 1]
[ts,ys] = ode45(f1,[0,200],[y10;y20]);
[ts1,ys1] = ode45(f2,[200,350],[y10;y20]);
[ts2,ys2] = ode45(f1,[350,500],[y10;y20]);
[ts3,ys3] = ode45(f2,[500,700],[y10;y20]);
[ts4,ys4] = ode45(f1,[700,900],[y10;y20]);
[ts5,ys5] = ode45(f1,[900,1000],[y10;y20]);
figure(7);
plot(ys(:,1),ys(:,2))
hold on;
plot(ys(1,1),ys(1,2),'bo') % starting point
hold on;
plot(ys(end,1),ys(end,2),'ks') % ending point
hold on;
plot(ys1(:,1),ys1(:,2))
hold on;
ylim([-10 100])
plot(ys1(1,1),ys1(1,2),'bo') % starting point
hold on;
plot(ys1(end,1),ys1(end,2),'ks') % ending point
hold on;
plot(ys2(:,1),ys2(:,2))
hold on;
plot(ys3(:,1),ys3(:,2))
hold on;
plot(ys4(:,1),ys4(:,2))
hold on;
plot(ys5(:,1),ys5(:,2))
xlabel('x1')
ylabel('x2')
legend('Phase Portrate x1,x2')
end
end
0 commentaires
Réponses (0)
Voir également
Catégories
En savoir plus sur Array Geometries and Analysis dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!