我有两条曲线(非线性),比如 C1 和 C2。彼此相交(可以不止一次)。我只需要使用 Scilab 找到这些交点。
问问题
2614 次
2 回答
1
在Open Source Engineering的 Scilab 教程中有两个例子:
在 Scilab 中求解非线性系统[PDF]。
于 2012-11-28T12:43:10.633 回答
0
你好 Scilab 和 Matlab 社区,
虽然这是一个古老的线程,但我会发布一个答案。由插值点给出的两条曲线的交点是一个标准问题。我必须为特定的应用程序解决它,在网络上发现了许多有用的灌输,但没有可行的解决方案。下面是函数 X_Crossing 应用于为该线程发明的非平凡情况(利萨如线和由 NaN 分隔的两个椭圆之间的一堆交点。)为方便起见,复制并粘贴到 Scilab 控制台中应该足以显示结果。
祝你好运
蔷薇花
对不起,我不明白你的代码块说明。抱歉,我没有找到预览按钮。
clear
clc()
mode(0)
format(10)
function [xC,yC,nC,dxC,dyC,duC,dvC]=X_Crossing(x,y,u,v,d_near,info)
//Alle Schnittpunkte zweier Polygonzüge (offen oder geschlossen))
//Selbstüberschneidungen sind keine Schnittpunkte!
//INPUT
//x.y. u,v: Koordinaten der beiden Kurvwen, 1 x nxy bzw. 1 x nuv
//Wenn möglich, alle Polygonkanten etwa gleich lang.
//d_near: Schnittpunkte werden nur in den Bereichen der Kurven
// gesucht, in denen die Kurven sich näher als d_near sind.
// d_near zu klein: Schnittpunkte werden übersehen.
// d_near unnötig groß: Rechnung dauert unnötig lange.
// Erster Versuch: d_near ca. dreifacher Punktabstand
//OUTPUT
//xC,yC: Koordinaten der Scnittpunkte, Zeilenvektoren
//nC: Anzahl der gefundenen Schnittpunkte
//info: 1: Hinweise, 0: keine Hinweise
//METHODE
// Scnittpunkte der Kantengeraden
warning('off')
Meldung=[' ']//Initialisierung
dxC=[];
dyC=[];
duC=[];
dvC=[];
nC=0//Initialisierung Anzahl Schnittpunkte
nxy=length(x)
nuv=length(u)
//Abschnitte in denen die Kurven nahe verlaufen
I=[];
for j=1:nuv
//Abstand des j-ten W-Punktes von SFK
d=sqrt([u(j)-x].^2+[v(j)-y].^2);
i_near=find(d<d_near);
I=[I i_near];//Indizes einer Untermenge von x,y
end//for j=1:nuv
I=unique(gsort(I));
I=I';
nI=length(I)
diffIgt1=find(diff(I)>1)';
Ia=[I(1); I(diffIgt1(1:$)+1)];
Ie=[I(diffIgt1(1:$));I($)];
nNahbereich=length(Ie);
//printf('\n nI: %f \n',nI)
//printf('\n Ia: %f \n',Ia)
//printf('\nnNahbereich: %f \n',nNahbereich)
if nNahbereich==0 then
if info
infotext=['Der Kurvenabstand ist überall > d_near.';
'Wenn Schnittpunkte existieren, d_near vergrößern!'];
Meldung=[Meldung;infotext]
printf('\n\n%4.0f Schnittpunkte gefunden.\n',nC)
printf('%s\n',Meldung)
end//if info
xC=%nan;
yC=%nan;
return
end// if nNahbereich= == 0 then
xC=[];//Schittpunktkoordinaten
yC=[];//Schittpunktkoordinaten
for k=1:nNahbereich
xNah=x(Ia(k):Ie(k));
yNah=y(Ia(k):Ie(k));
for i=1:length(xNah)-1
for j=1:nuv-1
dx=xNah(i+1)-xNah(i);
dy=yNah(i+1)-yNah(i);
du=-u(j+1)+u(j);
dv=-v(j+1)+v(j);
dxy=[dx;dy];
duv=[du;dv];
if ~isnan(dx)&~isnan(dy)&~isnan(du)&~isnan(dv)..
& dx*dv~=dy*du// nicht parallel
M=[dxy duv];//2 x 2
t=M\[u(j)-xNah(i);v(j)-yNah(i)];
if t(1)>=0 & t(1)<1 & t(2)>=0 & t(2)<1
xyC=[xNah(i);yNah(i)]+t(1)*dxy;
xC=[xC xyC(1)];
yC=[yC xyC(2)];
dxC=[dxC dx]
dyC=[dyC dy]
duC=[duC du]
dvC=[dvC dv]
end//if t(1)>=0 & ...
end//if ~isnan(dx)& ...
end//for j=1:nuv-1
end//for i=1:nxy-1
if length(xC)>0;
nC=length(xC)
// browsevar()
// abort
if info
printf('\n %5.0f Schnittpunkte gefunden',nC)
end//if info
end//if length(xC)>0;
if isempty(xC)|isempty(yC)
if info
infotext=['Wenn mehr existieren, d_near vergrößern!'];
Meldung=[Meldung; infotext]
printf('\n\n%4.0f Schnittpunkte gefunden.\n',nC)
printf('%s\n',Meldung)
end//if info
xC=%nan;
yC=%nan;
return
end//if isempty(xC)|isempty(yC)
end//for k=
endfunction//function [xC,yC,I,J]=X_Crossing(x,y,u,v)
if 1//just for stack overflow's forum parser
phiL=linspace(0,360,300);//Winkel für Lissajous-Figur
phiE=linspace(0,360,60);//Winkel für Ellipsen
//Zwei Ellipsen
x=-0.05+0.9*cosd(phiE);
x=[x %nan 0.5*x+0.5];
y=0.3+0.6*sind(phiE);
y=[y %nan 0.5*y-0.8];
u=cosd(3*phiL);
v=sind(5*phiL);
d_near=0.28;
info=1;
printf('\nBitte warten!\n')
printf(getversion()+' sucht Schnittpukte...')
xdel()
fig=figure('position',[450 -50 700 700],'background',8,'visible','off');
[xC,yC,nC,dxC,dyC,duC,dvC]=X_Crossing(u,v,x,y,d_near,info);
if ~isnan(xC)|~isnan(yC)
plot(xC,yC,'or')//Schnittpunktmarker
end//if ~isnan(xC)&~isnan(yC)
plot(x,y,u,v)//
// //Schnittpunkt-Tangenten
// plot([xC;xC+dxC], [yC;yC+dyC],'r','thickness',2)
// plot([xC;xC+duC], [yC;yC+dvC],'r','thickness',2)
ax=gca();
ax.isoview='on';
set(fig,'visible','on')
// browsevar()
end//if 1//just for stack overflow's forum parser
于 2020-03-10T20:27:04.033 回答