[总结]关于隐函数曲线和曲面画图的问题

xjzuo

已经多次看到论坛上有人问“隐函数曲线和曲面画图的问题”，特此作个小结,以便参考。
学习这些方法时,可以参考Matlab自带的help文件,以及论坛相关的讨论贴。
------这也正好是我在振动论坛的第1000贴,就算作为一个纪念吧。
%%%=======================================================================%%%
%%% 隐函数曲线画图的问题
%%% 方法有四：
1. 用ezplot
%%% 示例 画出隐函数z=x^2+y^2-5*x*y+1图形
%%%------------------------------------------%%%
ezplot('x^3+y^3-5*x*y+1',[-5,5])
%%%-----------------------------------------%%%

%%% 注：这里举个简单例子,再复杂的情形也只需套用即可
%%% 不过我还是希望读者自己也能多加思考,理解并创新

2. 用contour;
%%% 示例 用contour指令画出隐函数z=x^2+y^2-5*x*y+1图形
%%%------------------------------------------------%%%
clear all
x=linspace(-5,5);
y=linspace(-5,5);
[x1,y1]=meshgrid(x,y);
z1=x1.^3+y1.^3-5*x1.*y1+1;
contour(x1,y1,z1,[0 0],'b');
% axis image
%%%------------------------------------------------%%%

3. 直接利用Mathematica,Maple的函数ImplicitPlot函数画图，具体用法可参考其help文件中的例子。
%%% 示例
%%%------------------------------------------------%%%
<< Graphics`ImplicitPlot`
ImplicitPlot[x^3+y^3-5*x*y+1==0, {x, -5, 5}]
%%%------------------------------------------------%%%

4. 参考Mathematica, Maple，自己写一个函数实现。
%%% 这里贴一个mathworks网站上，外国人写的函数implot.m
%%% 其实也是借助了contour函数
%%% example
%%%------------------------------------------------%%%
clear all
f=inline('x^3+y^3-5*x*y+1','x','y');
implot(f,[-5 5 -5 5])
%%%
function implot(fun,rangexy,ngrid)
% Implicit plot function
% function implot(fun,rangexy,ngrid)
% fun is 'inline' function f(x,y)=0 (Note function written as equal to zero)
% rangexy =[xmin,xmax,ymin,ymax] range over which x and y is ploted default(-2*pi,2*pi)
% ngrid is the number of grid points used to calculate the plot,
% Start with course grid (ngrid =20) and then use finer grid if necessary
% default ngrid=50
%
% Example
% Plot y^3+exp(y)-tanh(x)=0
%
% write function f as an 'inline' function of x and y-- right hand side
% equal to zero
%
% f=inline('y^3+exp(y)-tanh(x)','x','y')
% implot(f,[-3 3 -2 1])
% A.Jutan UWO 2-2-98  ajutan@julian.uwo.ca
   if nargin == 1; % grid value and ranges not specified calculate default rangexy=[-2*pi,2*pi,-2*pi,2*pi];
   ngrid=50;
end
  if nargin == 2;  % grid value not specified
   ngrid=50;
end
  % get 2-D grid for x and y
  xm=linspace(rangexy(1),rangexy(2),ngrid);
ym=linspace(rangexy(3),rangexy(4),ngrid);
[x,y]=meshgrid(xm,ym);
fvector=vectorize(fun);% vectorize the inline function to handle vectors of x y
fvalues=feval(fvector,x,y); %calculate with feval-this works if fvector is an m file too
%fvalues=fvector(x,y); % can also calculate directly from the vectorized inline function
contour(x,y,fvalues,[0,0],'b-');% plot single contour at f(x,y)=0, blue lines
xlabel('x');ylabel('y');
grid
%%%--------------------------------------------------------------------------------------%%%

至于耦合方程的情形，可以利用Mathematica的函数ImplicitPlot实现,参考其help文件中的例子即明白。
前不久好像也讨论过该问题，可以参考之。

%%%=======================================================================%%%
%%% 隐函数曲面画图的问题
%%% 方法有四：
1. Maple有一个现成的函数implicitplot3d可用。
2. Mathematica调用程序包时，也可以利用函数ImplicitPlot3D实现隐函数曲面的绘制。
3. Mathematica中ContourPlot3D也可以画出隐函数曲面图,只是图形看上去不是那么漂亮。
4.目前Matlab并没有现成的函数可用,ezplot3只是画(参数)曲线图。但我们可以依循“利用等值线函数画隐函数曲线图形”的思路，利用现成的等值面函数来实现隐函数曲面的绘制。其实也可以看作是Mathematica中ContourPlot3D的类比。
%%% 举一个例子，原问题见：http://forum.vibunion.com/thread-54410-1-1.html
%%%---------------------------------------------------------------%%%
clear all
I0=10000;
sigma1 = 1;
sigma2 = 1;
MU = [0 0];
x=-10:0.1:10;y=-10:0.1:10;z=0:0.1:10;
[X,Y,Z]=meshgrid(x,y,z);
I1 =I0*1/(sqrt(2*pi)*sigma1*sigma2)*exp(-((Y+Z-MU(1)).^2/(4*sigma1^2)+(Z-Y-MU(2)).^2/(4*sigma2^2)...
    +(sigma1^2+sigma2^2)*X.^2/(2*sigma1^2*sigma2^2)));
p=patch(isosurface(X,Y,Z,I1,1));
set(p, 'FaceColor', 'red', 'EdgeColor', 'none');
daspect([1 1 1])
view(3)
camlight; lighting phong
%%%--------------------------------------------------------------%%%
%%%===============================================================%%%
当然，还有一些其它方法也可以实现以上问题的求解，例如神经网络等，
有兴趣的可以查找一下这方面的资料，这里就不一一列举了。

donkeyxu

原问题"[X,Y,Z] = meshgrid(-2:0.1:2,-2:0.1:2,-2:0.1:2);"

被改成x=-10:0.1:10;y=-10:0.1:10;z=0:0.1:10; [X,Y,Z]=meshgrid(x,y,z); 了啊

=========================
注: 请抓住本贴的目的，而不要本末倒置。
By xjzuo
=========================

[ 本帖最后由 xjzuo 于 2007-11-10 16:06 编辑 ]

donkeyxu

原帖由 donkeyxu 于 2007-11-10 14:54 发表

                               
登录/注册后可看大图

原问题"[X,Y,Z] = meshgrid(-2:0.1:2,-2:0.1:2,-2:0.1:2);"

被改成x=-10:0.1:10;y=-10:0.1:10;z=0:0.1:10; [X,Y,Z]=meshgrid(x,y,z); 了啊

=========================
注: 请抓住本贴的目的，而不要本末 ...

我是说这样改就能解决原来的问题了吗？我系统被破坏，matlab装不上，没法试。