forcal 发表于 2011-1-13 15:05:48

关于多值隐函数绘图问题

有关于多值隐函数绘图的问题向大家请教！

以这个函数为例（x取-8～8，y取-8～8）：f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;

可以求得x为任意值时对应的多个y值。绘图时，将每一组值(x,y)绘制为一个点，点数多了就会形成曲线，如下图：

但如何将相邻的两个点连成线呢？就是说如何判断两个点是相邻的呢？

例如，以下代码求得了x＝1.1及x＝1.2时对应的y值；x＝1.1的y值有6个，x＝1.2的y值也有6个，如何判断哪些点该连成线呢？!using("fcopt");
f(x,y)=(x^2+y^2)^3-36*(x^2-y^2)^2;
isolve;
isolve;1.1                     -5.413569621185505      0.
1.1                     -1.52029233271966         0.
1.1                     -0.8671650244626143       8.881784197001252e-016
1.1                     0.8671650244626143      8.881784197001252e-016
1.1                     1.52029233271966          0.
1.1                     5.413569621185504         0.
6.
1.2                     -5.274704716901752      0.
1.2                     -1.731609595222745      0.
1.2                     -0.9268324331975657       0.
1.2                     0.9268324331975657      0.
1.2                     1.731609595222745         0.
1.2                     5.274704716901752         0.
6.

wayne 发表于 2011-1-13 16:40:58

1# forcal
点画多了,越稠密,就越连续了,
================

或者求出一定数目的点,然后线性插值,样条插值,...............

wayne 发表于 2011-1-13 17:08:46

不知道你具体是在问什么

是关于画图引擎的渲染技巧吗
还是就题论题,只是画出该函数的 点图,然后标记一条直线?
===================================

前者,你可以参考gnuplot
如果是后者,你可以把方程化成极坐标的形式,就清晰了

forcal 发表于 2011-1-13 17:29:49

不知道你具体是在问什么

是关于画图引擎的渲染技巧吗
还是就题论题,只是画出该函数的 点图,然后标记一条直线?
===================================

前者,你可以参考gnuplot
如果是后者,你可以把方程化成极 ...
wayne 发表于 2011-1-13 17:08 http://bbs.emath.ac.cn/images/common/back.gif
画出函数的点图是比较容易的，我想把点连成一条折线（曲线）。
例如matlab等画这个图的曲线是连续的，但并不是点图。

＝＝＝＝＝＝＝＝＝
有些隐函数方程化成极坐标的形式是比较困难的。

wayne 发表于 2011-1-13 18:56:08

本帖最后由 wayne 于 2011-1-13 18:58 编辑

4# forcal
例如matlab等画这个图的曲线是连续的，但并不是点图。
:)
那你做一个测试,试试:
先产生一个点序列，100个 ，,比如   {x,sin(x)}
然后shuffle一下，把这100对的顺序全部搞乱，越乱越好，再plot一下，看看，..................

wayne 发表于 2011-1-13 19:10:21

4# forcal
如果你说的是contour plot,
那个实现起来要复杂一些, 你用google搜搜 contour plot Implementation

三天后我会回来的

wayne 发表于 2011-1-13 20:24:28

4# forcal
算了,再给你补一贴吧:

用OpenGL 画你给的那个函数:


代码:/* generates contours using marching squares */

/* region size */

#define X_MAX 8.0
#define Y_MAX 8.0
#define X_MIN -8.0
#define Y_MIN -8.0

/* number of cells */

#define N_X 500
#define N_Y 500

/* contour value */

#define THRESHOLD 0.0

#include <GL/glut.h>

void display()
{
   double f(double,double);
   int cell(double, double, double, double);
   void lines(int, int, int, double, double, double, double);

   double data;
   int i,j;
   int c;

   glClear(GL_COLOR_BUFFER_BIT);
/* form data array from function */

    for(i=0;i<N_X;i++) for (j=0;j<N_Y;j++)
       data=f(X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0), Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0));

    /* process each cell */

    for(i=0;i<N_X;i++) for (j=0;j<N_Y;j++)
    {
       c=cell(data, data, data, data);
       lines(c,i,j,data, data, data, data);
    }
    glFlush();
}

/* define function f(x,y)      */

double f(double x, double y)
{

   double a=36;

   /* Ovals of Cassini*/

   return (x*x+y*y)*(x*x+y*y)*(x*x+y*y)-a*(x*x-y*y)*(x*x-y*y);
}


/* define cell vertices */

int cell(double a, double b, double c , double d)
{
    int n=0;
    if(a>THRESHOLD) n+=1;
    if(b>THRESHOLD) n+=8;
    if(c>THRESHOLD) n+=4;
    if(d>THRESHOLD) n+=2;
    return n;
}


/* draw line segments for each case */

void lines(int num, int i, int j, double a, double b, double c, double d)
{
   void draw_one(int, int, int, double, double, double, double);
   void draw_adjacent(int, int, int, double, double, double, double);
   void draw_opposite(int, int, int, double, double, double, double);
   switch(num)
   {
   case 1: case 2: case 4: case 7: case 8: case 11: case 13: case 14:
      draw_one(num, i,j,a,b,c,d);
      break;
   case 3: case 6: case 9: case 12:
      draw_adjacent(num,i,j,a,b,c,d);
      break;
   case 5: case 10:
      draw_opposite(num, i,j,a,b,c,d);
      break;
   case 0: case 15: break;
   }
}

void draw_one(int num, int i, int j, double a, double b, double c, double d)
{
double x1, y1, x2, y2;
double ox, oy;
double dx, dy;
dx=(X_MAX-(X_MIN))/(N_X-1.0);
dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
    switch(num)
    {
    case 1: case 14:
      x1=ox;
      y1=oy+dy*(THRESHOLD-a)/(d-a);
      x2=ox+dx*(THRESHOLD-a)/(b-a);
      y2=oy;
      break;
    case 2: case 13:
      x1=ox;
      y1=oy+dy*(THRESHOLD-a)/(d-a);
      x2=ox+dx*(THRESHOLD-d)/(c-d);
      y2=oy+dy;
      break;
    case 4: case 11:
      x1=ox+dx*(THRESHOLD-d)/(c-d);
      y1=oy+dy;
      x2=ox+dx;
      y2=oy+dy*(THRESHOLD-b)/(c-b);
      break;
    case 7: case 8:
      x1=ox+dx*(THRESHOLD-a)/(b-a);
      y1=oy;
      x2=ox+dx;
      y2=oy+dy*(THRESHOLD-b)/(c-b);
      break;
    }
glBegin(GL_LINES);
    glVertex2d(x1, y1);
    glVertex2d(x2, y2);
glEnd();
}

void draw_adjacent(int num, int i, int j, double a, double b,
double c, double d)
{
double x1, y1, x2, y2;
double ox, oy;
double dx, dy;
dx=(X_MAX-(X_MIN))/(N_X-1.0);
dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
switch(num)
{
    case 3: case 12:
      x1=ox+dx*(THRESHOLD-a)/(b-a);
      y1=oy;
      x2=ox+dx*(THRESHOLD-d)/(c-d);
      y2=oy+dy;
      break;
    case 6: case 9:
      x1=ox;
      y1=oy+dy*(THRESHOLD-a)/(d-a);
      x2=ox+dx;
      y2=oy+dy*(THRESHOLD-b)/(c-b);
      break;
}
glBegin(GL_LINES);
    glVertex2d(x1, y1);
    glVertex2d(x2, y2);
glEnd();

}

void draw_opposite(int num, int i, int j, double a, double b,
double c, double d)
{
double x1,y1,x2,y2,x3,y3,x4,y4;
double ox, oy;
double dx, dy;
dx=(X_MAX-(X_MIN))/(N_X-1.0);
dy=(Y_MAX-(Y_MIN))/(N_Y-1.0);
ox=X_MIN+i*(X_MAX-(X_MIN))/(N_X-1.0);
oy=Y_MIN+j*(Y_MAX-(Y_MIN))/(N_Y-1.0);
switch(num)
{
    case5:
      x1=ox;
      y1=oy+dy*(THRESHOLD-a)/(d-a);
      x2=ox+dx*(THRESHOLD-a)/(b-a);
      y2=oy;
      x3=ox+dx*(THRESHOLD-d)/(c-d);
      y3=oy+dy;
      x4=ox+dx;
      y4=oy+dy*(THRESHOLD-b)/(c-b);
      break;
    case 10:
      x1=ox;
      y1=oy+dy*(THRESHOLD-a)/(d-a);
      x2=ox+dx*(THRESHOLD-d)/(c-d);
      y2=oy+dy;
      x3=ox+dy*(THRESHOLD-a)/(b-a);
      y3=oy;
      x4=ox+dx;
      y4=oy+dy*(THRESHOLD-b)/(c-b);
      break;
}
glBegin(GL_LINES);
    glVertex2d(x1, y1);
    glVertex2d(x2, y2);
    glVertex2d(x3, y3);
    glVertex2d(x4, y4);
glEnd();
}

void myReshape(int w, int h)
{
    glViewport(0, 0, w, h);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    if (w <= h) gluOrtho2D(X_MIN, X_MAX, Y_MIN * (GLfloat) h / (GLfloat) w, Y_MAX* (GLfloat) h / (GLfloat) w);
    else gluOrtho2D(X_MIN * (GLfloat) w / (GLfloat) h, X_MAX * (GLfloat) w / (GLfloat) h, Y_MIN, Y_MAX);
    glMatrixMode(GL_MODELVIEW);
}

int main(int argc, char **argv)
{
   glutInit(&argc, argv);
   glutInitWindowSize(500, 500);
   glutCreateWindow("contour plot");
   glutReshapeFunc(myReshape);
   glutDisplayFunc(display);
   glClearColor(0.0,0.0,0.0,1.0);
   glColor3f(1.0,1.0,1.0);
   glutMainLoop();
}

forcal 发表于 2011-1-13 20:24:43

4# forcal

:)
那你做一个测试,试试:
先产生一个点序列，100个 ，,比如   {x,sin(x)}
然后shuffle一下，把这100对的顺序全部搞乱，越乱越好，再plot一下，看看，..................
wayne 发表于 2011-1-13 18:56 http://bbs.emath.ac.cn/images/common/back.gif
我的意思是MATLAB提供了一个ezplot函数绘制隐函数图形，似乎不是用描点方式绘图的，它是如何绘制的？

对于单值函数f(x)=... ...；绘制曲线图是很容易的，将所有点按x的次序连起来即可。

但对于隐函数f(x,y)=0，一个x对应多个y值，该如何绘制曲线图？

forcal 发表于 2011-1-13 20:26:10

4# forcal
如果你说的是contour plot,
那个实现起来要复杂一些, 你用google搜搜 contour plot Implementation

三天后我会回来的
wayne 发表于 2011-1-13 19:10 http://bbs.emath.ac.cn/images/common/back.gif
英文不是很好，等你回来，呵呵。

forcal 发表于 2011-1-13 20:54:05

4# forcal
算了,再给你补一贴吧:

用OpenGL 画你给的那个函数:
wayne 发表于 2011-1-13 20:24 http://bbs.emath.ac.cn/images/common/back.gif
这个绘图原理与我所说的不一样，我原想先解出所有要绘制的点，然后连成线，看来比较复杂。

这个OpenGL 绘隐函数图原理能大致说一下吗？看不太懂。
页: [1] 2
查看完整版本: 关于多值隐函数绘图问题