11.matlab找两条离散曲线的交点
(2015-01-24 21:03:55)
标签:
股票 |
分类: Matlab |
引言
———————————————————分割线————————————————————————
对于1,我想说的是,如果你想要求得两曲线的精确交点,并且一个不漏,那就直接求解方程组,不用看本帖下文;
对于2,直线在Matlab里面是两个点确定,因此交点如果是一段线(无穷个点)的情况,可能只是显示两端点为交点;
对于3,很简单的例子,参数方程 x=cos(t),y=sin(t) (0<=t<=2*pi) 在数学分析(即连续空间)层面上是个圆,但是如果你在离散t的时候,间距比较大,那么最后Matlab绘制的图像不是圆,而是正多边形了。因此,此时我们讨论曲线交点是这个离散点连线的图形与其他图形的交点,而非圆与其他交点。这也是我在标题中加了“离散点连成”的修饰词,防止被误会。
对于4,既然是求曲线交点,那么本方法可以作为求方程组的近似解。当然,如果离散点够多,解的精确度可以保证,不过不能保证一个不漏。另外就是,对于一组离散点构成的曲线,很难知道它们的解析表达式,因此想通过非线性方程组求解的方法来求交点,就不大可能了(不过你可以用曲线拟合出函数解析式),因此,本帖的方法将会是一个较为有效求交点的方法。
下面将阐述我的方法以及给出例子代码。
下面是曲线y=cos(2*x).*exp(sin(x))与y2=sin(x).^2+cos(x)在[0:pi/18:2*pi]区间内的交点的代码:
注意:我没有写成接口的形式,虽然对于比那些较懒的人来说不太方便,但是这样做是为了让你能更好弄懂原理,并能自己改造代码。因此,下面的代码可以稍作修改,就能解决别的曲线求交点。这样,不愿思考的懒人就没法达到自己的目的了~
% 绘制两离散曲线的交点
% 注意:
%
%
%
%
%
%
clear;
debug=false; %关闭显示求交点过程
% 曲线1
x=0:pi/18:2*pi;
y=cos(2*x).*exp(sin(x));
% 曲线2
[x1 N]=sort(x);
% 下面几句代码在本个案下没有什么特殊作用,但是当出现参数方程的时候,下面的方法改动一下就会有用。
y1=sin(x1).^2+cos(x1); %用于作图
x2=x;
y2=sin(x).^2+cos(x); %用于寻点
h=plot(x1,y1,'b',x,y,'c');
y(abs(y)<=eps)=0; y2(abs(y2)<=eps)=0;%对于三角函数关于零点的部分处理,但是发现sin(k*pi)不一定全在eps范围内
cy=y2-y; %作差
%符号记录
pos=cy>0;
neg=cy<=0;
%确定变号位置
fro=diff(pos)~=0; %变号的前导位置
rel=diff(neg)~=0; %变号的尾巴位置
zpf=find(fro==1); %记录索引
zpr=find(rel==1)+1; %记录索引
zpfr=[zpf; zpr];
hold on
% 观看求交点过程
if debug, hp=plot(x(zpfr),y(zpfr),'r.-',x2(zpfr),y2(zpfr),'g.-'); end
%线性求交
x0=(x(zpr).*(y2(zpf)-y(zpf))-x(zpf).*(y2(zpr)-y(zpr)))./(y(zpr)+y2(zpf)-y(zpf)-y2(zpr));
y0=y(zpf)+(x0-x(zpf)).*(y(zpr)-y(zpf))./(x(zpr)-x(zpf));
if any(isnan(y0)), y0=y2(zpf); end
%加入已经判断为零的位置
x0=[x(abs(cy)<=eps) x0].';
y0=[y(abs(cy)<=eps) y0].';
hc=plot(x0,y0,'k.'); %绘制交点
if debug, legend([h;hc;hp(1);hp(end)],'C1','C2','交点','微线段1','微线段2',0); end
legend([h;hc],'C1','C2','交点',0)
xlabel('x'), ylabel('y'), zlabel('z');
title('平面曲线交点')
axis equal
hold off
disp('交点坐标[x,y]为:')
disp(unique([x0,y0],'rows')) %排除重复的点
经测试十几种奇怪的曲线相交(包括参数方程形式的曲线),目前发现上述代码的方法有四种情况会出现遗漏一两个交点。(其实上面代码本意是求显式函数的曲线交点,或者未知表达式的离散点曲线的交点,并未针对参数方程,隐函数方程做优化,但是可以凑合着用用。)