请教matlab的曲面间交线作图问题

sihaiwuya · 发表于 2005-11-24 22:57

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

请问如何用matlab画出两个曲面之间的交线。
例如，x^2+y^2+z^2=1与x^2-x+y^2=0的交线，最好能将这两个曲面画出来，交线用深颜色表示。

[ 本帖最后由 ChaChing 于 2009-2-27 17:40 编辑 ]

AaronSpark · 发表于 2005-11-25 01:01

套用一下下面的程序 clf, a=-20;eps0=1;
[x,y]=meshgrid(-10:0.2:10); %生成平面网格
v=[-10 10 -10 10 -100 100]; %设定空间坐标系的范围
colormap(gray) %将当前的颜色设置为灰色
z1=(x.^2-2*y.^2)+eps; %计算马鞍面函数z1=z1(x,y)
z2=a*ones(size(x)); %计算平面 z2=z2(x,y)
r0=abs(z1-z2)<=eps0; 
 %计算一个和z1同维的函数r0,当abs(z1-z2)<=eps时r0 =1；当abs(z1-z2)>eps0时，r0 =0。
 %可用mesh(x,y,r0)语句观察它的图形，体会它的作用，该方法可以套用。
zz=r0.*z2;xx=r0.*x;yy=r0.*y; %计算截割的双曲线及其对应的坐标
subplot(2,2,2), %在第2图形窗口绘制双曲线
 h1=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'+'); 
 set(h1,'markersize',2),hold on,axis(v),grid on
subplot(2,2,1), %在第一图形窗口绘制马鞍面和平面 
 mesh(x,y,z1);grid,hold on;mesh(x,y,z2); 
 h2=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'.'); %画出二者的交线
 set(h2,'markersize',6),hold on,axis(v),
for i=1:5 %以下程序和上面是类似的，通过循环绘制一系列的平面去截割马鞍面
 a=70-i*30; %在这里改变截割平面
 z2=a*ones(size(x)); r0=abs(z1-z2)<=1; zz=r0.*z2;yy=r0.*y;xx=r0.*x;
 subplot(2,2,3),
 mesh(x,y,z1);grid,hold on;mesh(x,y,z2);hidden off
 h2=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'.'); axis(v),grid
 subplot(2,2,4),
 h4=plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'o');
 set(h4,'markersize',2),hold on,axis(v),grid on
end

bainhome · 发表于 2005-11-25 01:51

转自《工程计算可视化与MATLAB实现》尚涛等编著武汉大学出版社 ------------动力与控制论坛.paradiseboy 两<PARSE>直线相交</PARSE> function [X,Y]=pll(X1,Y1,X2,Y2) % <PARSE>直线相交</PARSE>求交点 A1=Y1(1)-Y1(2); B1=X1(2)-X1(1); C1=Y1(2)*X1(1)-Y1(1)*X1(2); A2=Y2(1)-Y2(2); B2=X2(2)-X2(1); C2=Y2(2)*X2(1)-Y2(1)*X2(2); D=det([A1,B1;A2,B2]); X=det([-C1 B1;-C2 B2])/D; Y=det([A1 -C1;A2,-C2])/D; 调用格式： x1=[1 5];y1=[1 5];x2=[1 5];y2=[5,1]; [x,y]=pll(x1,y1,x2,y2); plot(x1,y1,'r'); hold on plot(x2,y2,'b'); plot(x,y,'ko'); %直线与多条<PARSE>直线相交</PARSE> xi=[1 2 3 4 5];yi=[2 6 3 6 1]; plot(xi,yi);hold on x1=[1 5];y1=[4 5];line(x1,y1); x=zeros(size(xi)); y=x; for i=1:5-1 x2=xi([i i+1]);y2=yi([i i+1]); [x,y]=pll(x1,y1,x2,y2); plot(x,y,'ro') end %直线与曲线相交 x=-8:0.1:8;y=x;[X,Y]=meshgrid(x,y); R=sqrt(X.^2+Y.^2)+eps;Z=sin(R)./R; contour(Z,3);hold on c=contour(Z,3); x=[0 360];y=[0 400]; y=(y(2)-y(1))/(x(2)-x(1))*(x-x(1))+y(1);z=[0 0]; line(x,y,z);c=c'; X=c(:,1);Y=c(:,2); r0=abs(Y-(y(2)-y(1))/(x(2)-x(1))*(X-x(1))+y(1))<=.93; zz=0;yy=r0.*Y;xx=r0.*X; plot(xx(r0~=0),yy(r0~=0),'r') %曲线与曲线相交 x=0:pi/400:2*pi; x=x'; y1=sin(pi*x);y2=cos(pi*x);plot(x,y1,x,y2);hold on r0=abs(y2-sin(pi*x))<=0.02; yy=r0.*y1;xx=r0.*x;plot(xx(r0~=0),yy(r0~=0),'r.') 直线与曲面相交 x=-8:0.3:8;y=x;[X,Y]=meshgrid(x,y); Z=X.^2+Y.^2; mesh(X,Y,Z);hold on x=[-10 10];y=[-10 3];z=[30 35];line(x,y,z); r0=(abs(Y-y(1)-(y(2)-y(1))/(x(2)-x(1))*(X-x(1)))<=0.45)&... (abs(Z-z(1)-(z(2)-z(1))/(x(2)-x(1))*(X-x(1)))<0.45)&... (abs(Y-y(1)-(y(2)-y(1))/(z(2)-z(1))*(Z-z(1)))<=0.45); zz=r0.*Z;yy=r0.*Y;xx=r0.*X; plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'r*') 平面与曲面相交 x=-8:0.1:8; y=x; [X,Y]=meshgrid(x,y); Z1=2*ones(size(X)); Z2=X.^2-Y.^2; mesh(X,Y,Z1); hold on mesh(X,Y,Z2); r0=(abs(Z1-Z2)<=.65); zz=r0.*Z1;yy=r0.*Y;xx=r0.*X; plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'k*') clc disp('观察曲面后，按任意键画交线'); pause clf plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'k*'); %曲面与多个截平面相交 y=-10:0.5:10; z=y; [Z,Y]=meshgrid(z,y); X=Z; X1=0*ones(size(Z)); X2=3*ones(size(Z)); X3=-3*ones(size(Z)); Z4=(X.^2-Y.^2)/10; mesh(X1,Y,Z);hold on mesh(X2,Y,Z) mesh(X3,Y,Z); mesh(X,Y,Z4); r1=(abs(X1-X)<0.05); r2=(abs(X2-X)<0.05); r3=(abs(X3-X)<0.05); zz1=r1.*Z4;yy1=r1.*Y;xx1=r1.*X; zz2=r2.*Z4;yy2=r1.*Y;xx2=r1.*X; zz3=r3.*Z4;yy3=r1.*Y;xx3=r1.*X; plot3(xx1(r1~=0),yy1(r1~=0),zz1(r1~=0),'k*'); plot3(xx2(r2~=0),yy2(r2~=0),zz2(r2~=0),'k*'); plot3(xx3(r3~=0),yy3(r3~=0),zz3(r3~=0),'k*'); colormap(hsv) clc; disp('观察曲面后，按任意键画交线'); hold off 平面与曲面相交 y=-8:0.4:8; z=y; [Z,Y]=meshgrid(z,y); X=Z; X1=zeros(size(Z)); Z2=zeros(size(Z)); Z3=(X.^2-Y.^2)/10; mesh(X1,Y,Z); hold on mesh(X,Y,Z2); mesh(X,Y,Z3); r1=(abs(X1-X)<0.05); r2=(abs(Z3-Z2)<0.05); r3=(abs(X1-X)<0.05)&(abs(Z-Z2)<=0.05); zz1=r1.*Z3;yy1=r1.*Y;xx1=r1.*X; zz2=r2.*Z3;yy2=r2.*Y;xx2=r2.*X; zz3=r3.*Z;yy1=r3.*Y;xx1=r3.*X1; plot3(xx1(r1~=0),yy1(r1~=0),zz1(r1~=0),'k*'); plot3(xx2(r2~=0),yy2(r2~=0),zz2(r2~=0),'k*'); plot3(xx3(r3~=0),yy3(r3~=0),zz3(r3~=0),'k*'); colormap(hsv);

sihaiwuya · 发表于 2005-11-25 17:59

谢谢了.
能给出我问的问题的解法吗？比如说，如何画出柱面x^2+y^2=1的图形。

[ 本帖最后由 ChaChing 于 2009-2-27 17:37 编辑 ]

tongjie · 发表于 2005-12-5 18:47

不错的方法

fifor · 发表于 2009-2-27 14:02

这个问题学习了，很受益

friendchj · 发表于 2009-2-27 15:01

参考：http://forum.vibunion.com/forum/ ... p;extra=&page=1

friendchj · 发表于 2009-2-27 15:23

clear
clc
[X,Y,Z]=sphere;
[X1,Y1,Z1]=cylinder;
X1=.5*X1+.5;
Y1=.5*Y1;
surf(X1,Y1,Z1*4-2)
hold on
surf(X,Y,Z)
x1=0:0.01:1;
y1=sqrt(x1-x1.^2);
y2=-y1;
z1=sqrt(1-x1);
z2=-z1;
plot3(x1,y1,z1,'r','LineWidth',5)
plot3(x1,y1,z2,'r','LineWidth',5)
plot3(x1,y2,z1,'r','LineWidth',5)
plot3(x1,y2,z2,'r','LineWidth',5)

复制代码

账号		自动登录	找回密码
密码			我要加入

[工具箱] 请教matlab的曲面间交线作图问题

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

3x

评分