آموزش نرم‌افزار متلب (MATLAB)تصویر فاز (Phase Portrait) مجموعه‌ای از تراژکتوری‌ها (Trajectory) هستند که معمولا برای سیستم‌های دو بعدی و یا سیستم‌هایی که شامل دو متغییر حالت هستند رسم می‌شود. به کمک Phase Portrait می‌توان به سادگی نقاط تعادل سیستم را در ناحیه رسم شده مشاهده کرد و پایداری یا ناپایداری آنها را تشخیص داد. رسم Phase Portrait در علوم مختلف از جمله کنترل اهمیت ویژه‌ای دارد. با رسم مجموعه‌ای از تراژکتوری‌ها با شرایط اولیه مختلف می‌توان به تصویر فاز سیستم دست یافت که در ادامه نحوه رسم Phase Portrait در نرم‌افزار متلب شرح داده خواهد شد.
فهرست مطالب
رسم Phase Portrait در نرم‌افزار متلب (MATLAB)
در این روش برای رسم Phase Portrait به سه تابع نیاز داریم. تابع اول برای معرفی سیستم مورد استفاده قرار می‌گیرید. دو تابع دیگر برای رسم و تعیین مسیر تراژکتوری‌ها تعریف می‌شوند و در نهایت با استفاده از این توابع تصویر فاز سیستم رسم می‌شود.
تابع معرفی سیستم در متلب
ابتدا باید سیستم مورد نظر خود را در قالب یک تابع معرفی نمایید. توجه داشته باشید این آموزش تنها برای رسم تصویر فاز سیستم‌های دو بعدی و یا به عبارت دیگر سیستم‌هایی که تنها شامل دو متغییر حالت هستند، قابل استفاده است.
مثال اول: معادله فضای حالت سیستم غیرخطی به صورت زیر است (تنها شامل دو متغییر حالت):
System\,\,\text{1\,\,:\,\,}\dot{x}_1=2x_1-x_1x_2\,\,,\,\,\dot{x}_2=2x_{1}^{2}-x_2
تابع این سیستم در متلب به صورت کد زیر تعریف می‌شود:
|
1 2 3 4 5 6 7 |
function [xdot] = system1(t,x) if nargin==1 x=t; end xdot(1,1)=2*x(1)-x(1)*x(2); xdot(2,1)=2*x(1)^2-x(2); end |
مثال دوم: معادله فضای حالت سیستم غیرخطی وندرپل (Van der Pol) به صورت زیر است:
Van\,\,der\,\,Pol\,\,\left( \mu =2 \right) \,\,:\,\,\dot{x}_1=x_2\,\,,\,\,\dot{x}_2=-x_1-2\left( x_{1}^{2}-1 \right) x_2
تابع سیستم وندرپل در متلب به صورت کد زیر تعریف می‌شود:
|
1 2 3 4 5 6 7 |
function [xdot] = vanderpol(t,x) if nargin==1 x=t; end xdot(1,1)=x(2); xdot(2,1)=-x(1)-2*(x(1)^2-1)*x(2); %for mu=2 end |
توابع arrowh و PhasePlane
برای رسم و تعیین مسیر تراژکتوری‌ها و نمایش Field Vector دو تابع با نام‌های arrowh و PhasePlane تعریف می‌شود که تنها کافیست آنها را با نام arrowh و PhasePlane ذخیره نمایید. تابع arrowh به صورت کد زیر تعریف می‌شود:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 |
function handle = arrowh(x,y,clr,ArSize,Where) %-- errors if nargin < 2 error('Please give enough coordinates !'); end if (length(x) < 2) || (length(y) < 2), error('X and Y vectors must each have "length" >= 2 !'); end if (x(1) == x(2)) && (y(1) == y(2)), error('Points superimposed - cannot determine direction !'); end if nargin <= 2 clr = 'b'; end if nargin <= 3 ArSize = [100,100]; end handle = []; %-- check if variables left empty, deal with ArSize and Color if isempty(clr) clr = 'b'; nonsolid = false; elseif ischar(clr) if strncmp('e',clr,1) % for non-solid arrow heads nonsolid = true; clr = clr(2); else nonsolid = false; end elseif isvector(clr) if length(clr) == 4 && clr(1) == 0 % for non-solid arrow heads nonsolid = true; clr = clr(2:end); else nonsolid = false; end else error('COLOR argument of wrong type (must be either char or vector)'); end if nargin <= 4 if (length(x) == length(y)) && (length(x) == 2) Where = 100; else Where = 50; end end if isempty(ArSize) ArSize = [100,100]; end if length(ArSize) == 2 ArWidth = 0.75*ArSize(2)/100; % .75 to make arrows slimmer else ArWidth = 0.75; end ArSize = ArSize(1); %-- determine and remember the hold status, toggle if necessary if ishold, WasHold = 1; else WasHold = 0; hold on; end %-- start for-loop in the case that several arrows are wanted for Loop = 1:length(Where), %-- vectors "longer" than 2 => we're dealing with time series if (length(x) == length(y)) && (length(x) > 2), j = floor(length(x)*Where(Loop)/100); %-- determine that location if j >= length(x), j = length(x) - 1; end if j == 0, j = 1; end x1 = x(j); x2 = x(j+1); y1 = y(j); y2 = y(j+1); else %-- just two points are given - take those x1 = x(1); x2 = (1-Where/100)*x(1)+Where/100*x(2); y1 = y(1); y2 = (1-Where/100)*y(1)+Where/100*y(2); end %-- get axe ranges and their norm OriginalAxis = axis; Xextend = abs(OriginalAxis(2)-OriginalAxis(1)); Yextend = abs(OriginalAxis(4)-OriginalAxis(3)); %-- determine angle for the rotation of the triangle if x2 == x1, %-- line is vertical, no need to calculate slope if y2 > y1, p = pi/2; else p= -pi/2; end else %-- line is not vertical, go ahead and calculate slope %-- using normed differences (seems better) m = ( (y2 - y1)/Yextend ) / ( (x2 - x1)/Xextend ); if x2 > x1, %-- now calculate the resulting angle p = atan(m); else p = atan(m) + pi; end end %-- the arrow is made of a transformed "template triangle". %-- it will be created, rotated, moved, resized and shifted. %-- the template triangle (it points "east", centered in (0,0)): xt = [1 -sin(pi/6) -sin(pi/6)]; yt = ArWidth*[0 cos(pi/6) -cos(pi/6)]; %-- rotate it by the angle determined above: xd = []; yd = []; for i=1:3 xd(i) = cos(p)*xt(i) - sin(p)*yt(i); yd(i) = sin(p)*xt(i) + cos(p)*yt(i); end %-- move the triangle, its "head" lays in (0,0): xd = xd - cos(p); yd = yd - sin(p); %-- stretch/deform the triangle to look good on the current axes: xd = xd*Xextend*ArSize/10000; yd = yd*Yextend*ArSize/10000; %-- move the triangle to the location where it's needed xd = xd + x2; yd = yd + y2; %-- draw the actual triangle handle(Loop) = patch(xd,yd,clr,'EdgeColor',clr); if nonsolid, set(handle(Loop),'facecolor','none'); end end % Loops %-- restore original axe ranges and hold status axis(OriginalAxis); if ~WasHold, hold off end end |
تابع PhasePlane نیز به صورت کد زیر تعریف می‌شود:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 |
function []=PhasePlane(system,tspan,icond,varargin) hy=0.5; hx=0.5; Ylim=[-5,5]; Xlim=[-5,5]; scale=1; PlotSingularPoints=false; ArrowHeads=false; ArrowSize=2; if ~isempty(varargin) for ii=1:length(varargin) if ischar(varargin{ii}) switch varargin{ii} case 'scale' scale=varargin{ii+1}; case 'Xlim' Xlim=varargin{ii+1}; case 'Ylim' Ylim=varargin{ii+1}; case 'hx' hx=varargin{ii+1}; case 'hy' hy=varargin{ii+1}; case 'PlotSingularPoints' % if isboolean(varargin{ii+1}) PlotSingularPoints=varargin{ii+1}; % else % error('the property PlotSingularPoints takes only 0 or 1 values'); % end case 'ArrowHeads' ArrowHeads=varargin{ii+1}; case 'ArrowSize' ArrowSize=varargin{ii+1}; otherwise warning(['Undefined property',varargin{ii}]); end end end end clear varargin; % solving system equations for all defined initial conditions for n=1:length(icond) [T,x]=ode23(system,tspan,icond{n}); clear T; x1=x(:,1); x2=x(:,2); hold on; plot(x1,x2,'k') if ArrowHeads index1=ceil(mean(0:length(x))); index2=ceil(mean(0:index1)); arrowh([x1(index1), x1(index1+1)], [x2(index1), x2(index1+1)],'r',ArrowSize*100); arrowh([x1(index2), x1(index2+1)], [x2(index2), x2(index2+1)],'r',ArrowSize*100); end end % quiver plot and singularities determination x_vec=Xlim(1):hx:Xlim(2); y_vec=Ylim(1):hy:Ylim(2); [X,Y]=meshgrid(x_vec,y_vec); U=zeros(size(X)); V=zeros(size(X)); for n=1:numel(X) Out=system([X(n),Y(n)]); U(n)=Out(1); V(n)=Out(2); if PlotSingularPoints if mod(n,10)==0 Singular=fsolve(system,[X(n),Y(n)]); if ~isnan(Singular) hold on; plot(Singular(1),Singular(2),'k*','LineWidth',1.2); end end end end hold on; quiver(X,Y,U,V,scale,'b') grid minor; box on; plot([x_vec(1),x_vec(end)],[0,0],'k-.',[0,0],[y_vec(1),y_vec(end)],'k-.'); xlim([x_vec(1),x_vec(end)]); ylim([y_vec(1),y_vec(end)]); xlabel('$x_1$','interpreter','latex','FontSize',14); ylabel('$x_2$','interpreter','latex','FontSize',14); end |
رسم تصویر فاز سیستم (Phase Portrait)
پس از ذخیره سه تابع معرفی شده، نوبت به رسم تصویر فاز سیستم (Phase Portrait) می‌رسد. هر یک از مثال‌های بالا توسط کد متفاوتی رسم می‌شود تا نحوه عملکرد توابع با ورودی‌های مختلف به خوبی نمایش داده شود. تابع PhasePlane سه ورودی ضروری دارد که عبارت است از سیستم، بازه زمانی و شرایط اولیه. سیستم قبلا به صورت تابع تعریف شده‌است که نام تابع پس از علامت @ قرار می‌گیرد. تراژکتوری‌ها در یک بازه زمانی رسم می‌شوند که این بازه زمانی با نام tspan مشخص شده‌است. بازه زمانی با توجه به نوع سیستم می‌تواند مقادیر مختلفی داشته باشد. شرایط اولیه نیز با متغییر icond وارد تابع PhasePlane می‌شود که در مسائل مختلف می‌تواند مقادیر مختلفی را به خود اختصاص دهد.
با استفاده از کد زیر می‌توان تصویر فاز سیستم مثال اول (System 1) را رسم کرد:
|
1 2 3 4 5 |
clc; clear; close all; tspan=[0,50]; icond={[0.01, -5], [-0.01,-5], [0.01,0], [-0.01,0]}; figure; PhasePlane(@system1,tspan,icond); |
تصویر فاز سیستم (Phase Portrait) برای مثال اول (System 1) به صورت شکل زیر رسم می‌شود:
همچنین با استفاده از کد زیر می‌توان تصویر فاز سیستم مثال دوم (Van der Pol) را رسم کرد:
|
1 2 3 4 5 |
clc; clear; close all; figure; tspan=[0,50]; icond={[0.1, 0.1], [5,5], [-5,5]}; PhasePlane(@vanderpol,tspan,icond,'scale',2,'ArrowSize',3,'ArrowHeads',true,'PlotSingularPoints',true); |
تصویر فاز سیستم وندرپل به صورت شکل زیر رسم می‌شود:
توجه داشته باشید که تراژکتوری‌ها را می‌توان از هر شرایط اولیه دلخواهی رسم کرد. از دستور حلقه می‌توان برای تولید شرایط اولیه و رسم تراژکتوری ها استفاده نمود.
دیدگاه بگذارید