В MATLAB для закрашивания области под кривой Y существует функция area(X,Y). Однако при построении графиков иногда требуется заштриховать (а не закрасить) область прямыми линиями, для чего в MATLAB нет встроенного инструмента. Для решения этой проблемы мной была написана функция plotShady(x1,y1,x2,y2,k0,r,color,lineW), которая покрывает область параллельными линиями между двумя кривыми: первая кривая - x1 и y1, вторая кривая - x2 и y2. Можно задавать их наклон (k0), толщину (lineW), цвет (color) и расстояние (r) между ними.
Например, чтобы покрыть область сеткой, необходимо вызвать функцию plotShady() два раза: сначала горизонтальные линии, т. е. с параметром k==0, затем вертикальные, т. е. с параметром k==NaN (рисунок 1).
x1=[0:10];x2=x1;
y1=sin(x1);y2=cos(x2);
plot(x1,y1);hold on;plot(x2,y2);
hold on;plotShady(x1,y1,x2,y2,0,0.1);
hold on;plotShady(x1,y1,x2,y2,NaN,0.7);
Следующий код заштриховывает область между двумя кривыми, у которых начало и конец не совпадают.
x1=[-5:5];x2=[-6:4];
y1=sin(x1)+5;y2=cos(x2)-2;
figure
plot(x1,y1);hold on;plot(x2,y2);
hold on;plotShady(x1,y1,x2,y2,-0.5,0.7);
Если нам нужно заштриховать область только под одной кривой, то значения второй кривой y2 умножаем на 0, как представлено в коде ниже.
x1=[-5:5];
y1=sin(x1)+4;
figure
plot(x1,y1,'-o','LineWidth',1.5);
hold on;plotShady(x1,y1,x1,y1*0,0.5,0.7);
Результат на рисунке 3.
Вот код функции plotShady():
function plotShady(x1,y1,x2,y2,k0,r,color,lineW)
% Отрисовка штриховых областей
% k=NaN - вертикальные прямые
% k=0 - горизонтальные прямые
% r – расстояние между линиями
if nargin<7 color=[127 127 127]/255; end
if nargin<8 lineW=0.5; end
%% =======корректировка исходных данных=======
if isnan(k0)==1
k=0;
else
k=k0;
end
x=sort(unique([x1,x2])); m1=0; m2=0;
if x1(1)>x2(1) y1=[y2(1),y1];x1=[x2(1),x1]; m1=1; end
if x1(1)<x2(1) y2=[y1(1),y2];x2=[x1(1),x2]; m1=1; end
if x1(end)<x2(end) y1=[y1,y2(end)];x1=[x1,x2(end)]; m2=1; end
if x1(end)>x2(end) y2=[y2,y1(end)];x2=[x2,x1(end)]; m2=1; end
y1=My_X_Y_X1(x1,y1,x,0);
y2=My_X_Y_X1(x2,y2,x,0);
y1_0=y1;y2_0=y2;
%% ====поиск точек пересечения исходных линий=======
dd=sign(y1-y2);
dd=[dd(1:end-1)-dd(2:end),0];
N=length(dd);
ind=[1:N];
ind=ind(dd~=0);
k1_0=(y1(ind+1)-y1(ind))./(x(ind+1)-x(ind));b1_0=y1(ind)-k1_0.*x(ind);
k2_0=(y2(ind+1)-y2(ind))./(x(ind+1)-x(ind));b2_0=y2(ind)-k2_0.*x(ind);
x_0=-(b2_0-b1_0)./(k2_0-k1_0);y_0=k1_0.*x_0+b1_0;
x=[x,x_0];y1_0=[y1_0,y_0];y2_0=[y2_0,y_0];
[x,ind]=sort(x);y1_0=y1_0(ind);y2_0=y2_0(ind);
y1=max([y1_0;y2_0]);y2=min([y1_0;y2_0]);
yyy=y1-y2;yyy(yyy<1e-9)=0;
%% ================
% figure;hold on;plot(x,y1,'-g*');hold on;plot(x,y2,'-r*');
%% =============================
dy=abs(max(y1)-min(y2));
% if k0~=0 && isnan(k0)==0 dx=dy/k; else dx=dy*k; end
if abs(k0)<1 && k0~=0 && isnan(k0)==0 dx=dy/k; else dx=dy*k; end
dx=abs(dx);
dn=round(dx/r);
y1=[y1(1),y1,y1(end)];
y2=[y2(1),y2,y2(end)];
x0=x;x=[x(1)-dx,x,x(end)+dx];
n0=length(x);
if k0~=0
n=abs(round((x(end)-x(1))/r))+2*dn;
else
n=abs(round((max(y1)-min(y2))/r))+2*dn;
end
if k0>=0 || isnan(k0)==1 y0=min(y2); end
if k0<0 y0=max(y1); end
b0=r*k;
if k0==0 b0=r; end
b=k*x(1)+y0;
b=b+[0:n-1]*b0;
%% ===============
k1=(y1(2:end)-y1(1:end-1))./(x(2:end)-x(1:end-1));b1=y1(1:end-1)-k1.*x(1:end-1);
k2=(y2(2:end)-y2(1:end-1))./(x(2:end)-x(1:end-1));b2=y2(1:end-1)-k2.*x(1:end-1);
kk1=repmat(k1,n,1);
kk2=repmat(k2,n,1);
l1=repmat(b1,n,1);
l2=repmat(b2,n,1);
bb=repmat(b',1,n0-1);
if isnan(k0)==1 % вертикальные прямые
x00=x0(1)+[0:n-1]*r;
X1=repmat(x00',1,n0-1);X2=X1;
Y1=kk1.*X1+l1;
Y2=kk2.*X2+l2;
elseif k0==0 % горизонтальные прямые
Y1=bb;X1=(Y1-l1)./kk1;
Y2=bb;X2=(Y2-l2)./kk2;
else
Y1=(l1*k-b'*k1)./(k-kk1);X1=(Y1-bb)/k;
Y2=(l2*k-b'*k2)./(k-kk2);X2=(Y2-bb)/k;
end
%% ==========================
ind1=zeros(size(X1));
ind2=zeros(size(X2));
for i=1:n0-1
ind1(X1>=x(i) & X1<=x(i+1) & kk1==k1(i) & l1==b1(i))=1;
ind2(X2>=x(i) & X2<=x(i+1) & kk2==k2(i) & l2==b2(i))=1;
end
Y1=Y1(ind1~=0);X1=X1(ind1~=0);bb1=bb(ind1~=0);
Y2=Y2(ind2~=0);X2=X2(ind2~=0);bb2=bb(ind2~=0);
Y=[Y1',Y2'];X=[X1',X2'];
% hold on;plot(X,Y,'mo')
%% ==================================
ind1=X1>=x0(1) & X1<=x0(end);Y1=Y1(ind1);X1=X1(ind1);bb1=bb1(ind1);
ind2=X2>=x0(1) & X2<=x0(end);Y2=Y2(ind2);X2=X2(ind2);bb2=bb2(ind2);
Y=[Y1',Y2'];X=[X1',X2'];bb=[bb1',bb2'];
XX=[];YY=[];
NN=[1:length(x0)+1];
if isnan(k0)==1 % вертикальные прямые
x0=x(1)+[0:n-1]*r;
for i=2:length(x0)
xx=[x0(i), x0(i)];
yy=Y(X==x0(i));
myy=max(abs(yy));
if myy~=0
yy=round(yy/myy*(1e+6))/(1e+6)*myy; % две близкие точки из-за погрешности
end
yy=unique(yy);
yy=[yy(1) yy(end)];
YY=[YY;yy];XX=[XX;xx];
end
else
for i=1:length(b)
xx=X(bb==b(i));
yy=Y(bb==b(i));
if m1==0
yy0=k*x0(1)+b(i);
if yy0<=y1(1) && yy0>=y2(1)
xx=[xx,x0(1)];yy=[yy,yy0];
end
end
if m2==0
yy0=k*x0(end)+b(i);
if yy0<=y1(end) && yy0>=y2(end)
xx=[xx,x0(end)];yy=[yy,yy0];
end
end
myy=max(abs(yy));
if myy~=0 yy=round(yy/myy*(1e+6))/(1e+6)*myy; end
mxx=max(abs(xx));
if mxx~=0 xx=round(xx/mxx*(1e+6))/(1e+6)*mxx; end
[xx,ind]=unique(xx);
yy=yy(ind);
% ===области, где функции совпадают
ax=xx;ay=yy;
for kk=1:length(xx)
n0=NN;
[xx0,ind]=sort([x0,ax(kk)]);
ty=[yyy,1];ty=ty(ind);
ind=n0(xx0==ax(kk));
if sum(ty(ind))==0
yy(xx==ax(kk))=[];
xx(xx==ax(kk))=[];
end
end
% ===
lxx=length(xx);
if mod(lxx,2)==0
xx=reshape(xx,[2,lxx/2])';
yy=reshape(yy,[2,lxx/2])';
XX=[XX;xx];
YY=[YY;yy];
end
end
end
% Y=[Y1',Y2'];X=[X1',X2'];
% hold on;plot(X,Y,'mo')
plot(XX',YY','Color',color,'LineWidth',lineW)
Функция My_X_Y_X1():
function [y1,k,b]=My_X_Y_X1(x,y,x1,sw)
%Функция «Линейная интерполяция»: подпрограмма для сопоставления исходных
%данных Y заданному массиву X1
[x,cc]=sort(x);y=y(cc);
if x(1)>x1(1) x(1)=x1(1); warning('Исходный x меньше текущего x1!'); end
if x(end)<x1(end) x(end)=x1(end); warning('Исходный x меньше текущего x1!'); end
% x1(1)=x(1);x1(end)=x(end);
k=(y(2:end)-y(1:end-1))./(x(2:end)-x(1:end-1));
b=y(1:end-1)-k.*x(1:end-1);
[x2,a]=sort([x1,x]);
if length(x1)>1
n=[1:length(x2)];
n0=zeros(1,length(x2));
n0(a>(length(x1)))=[1:length(x)];
ind=[];n00=n0;w=1;
while isempty(n0)==0
if length(n0)>1 && (n0(1)~=0 & n0(2)==0)
w=n0(1);
end
n0(1)=[]; ind=[ind,w];
end
ind(n00~=0)=[];
else
ind=find(x2==x1)-1;
if x1==x(1) ind=1; end
if x1==x(end) ind=length(x)-1; end
end
y1=k(ind).*x1+b(ind);
if sw~=0
figure;plot(x,y,'-r*',x1,y1,'--ok'); legend({'Old points' 'New points'});
end
Надеюсь, кому-то пригодится данное решение.
