如果用M代码提取穿越次数呢？

ppp546

如图所示，是一个二阶振荡曲线，在Workspace有两个一维向量，tout<4000*1>,yout<4000*1>,其中第399个值是0.9999，第400个值是1.0002，第401个值为1.003，我有什么办法能提取出来穿越1的次数呢？

iwbm

遍历 每两个相邻的数，看是否包含1？

xiezhh

我写一个例子，振荡曲线方程为：
y = 1-0.5*x.*cos(2*pi*x);
运行下面代码求穿过直线y=1的次数和交点

1. x=linspace(-8,0,4000);
   
2. y1 = 1-0.5*x.*cos(2*pi*x);
   
3. y2 = ones(size(y1));
   
4. y1_y2 = diff((y1-y2)>=0);
   
5. id = find(y1_y2);
   
6. cross_num = length(id)  %穿越次数
   
7. plot(x,y1,'r')
   
8. hold on
   
9. %求交点
   
10. xcross = [];
    
11. ycross = [];
    
12. if ~isempty(id)
    
13.     fenmu = y2(id)-y2(id+1)+y1(id+1)-y1(id);
    
14.     xcross = (x(id+1).*(y2(id)-y1(id))+x(id).*(y1(id+1)-y2(id+1)))./fenmu;
    
15.     ycross = (y1(id+1).*y2(id)-y1(id).*y2(id+1))./fenmu;
    
16. end
    
17. [xcross' ycross']
    
18. plot(xcross, ycross, 'o')

复制代码

效果图如下

ppp546

我也写了一个代码，可能运行效率没有楼上的高
for i=1:(length(y1)-1)
        if y1(i)<1
                if y1(i)(i+1)>=1
                        N=N+1;
                end;
        end;
end;
N=N/2;

ppp546

谢谢xiezhh

xiezhh

我觉得用matlab求曲线交点的时候，最正常的思路就是解方程，若求解析解，很可能出现无解，若求数值解，又跟初值的选取有很大关系。以前也有很多求曲线交点，线面交点的讨论，比如《工程计算可视化与MATLAB实现》（尚涛编著）里有很多现成的程序，那里的做法基本都是设定一个误差，让差的绝对值小于这个误差来求交点，对有些情况，结果很不精确。
        我的思路也是先将曲线数据离散化（变成等长的向量），然后确定两条曲线互穿的位置（即下标），就可以在两条曲线上分别找到相邻的两个点，从而确定两条直线，求这两条直线的交点即可得两条曲线的交点，运算都是基于向量化运算，避免了循环。当然这种作法求出的交点也会有一定误差，但是从图上看不出来，即使一直放大图像，我想原因大家应该能够明白。

ppp546

我明白xiezhh 的意思，通过等长离散化就可以求出曲线交战数值，而且也不需要编程。只要采样点足够多，就可保证交点解的准确性。
谢谢！