先把标题放在这里,内容我慢慢准备,一点点添加1.
一个简单的例子
从最简单的情况出发,连续可微的封闭曲线为边界的,最简单例子是”圆”.
——“什么？你忽悠我？！”
为了不至于激起不满，我临时决定把这个简单的例子更改为椭圆。当然，如果椭圆的长短轴跟直角坐标系的坐标轴平行，其实也还简单，所以让它再旋转一个角度 θ 。所以我选择了如下参数方程的椭圆：
{x=y=5costcosθ−3sintsinθ3cosθsint+5costsinθ
为了把曲线画出来, 让 θ=30o , 从而其参数方程变为:
⎧⎩⎨⎪⎪⎪⎪x=y=12(53√cost−3sint)12(5cost+33√sint)
ParametricPlot[RotationTransform[Pi/6][curve], {t, 0, 2 Pi}] /. Line[x_] :> {Red, Line[x]}
要判断任给点 P=(x0,y0) 是不是在图中的浅蓝色区域:
不知道常用的都是什么办法? 2.
.
至少对圆的情形,判断的方法是到圆心的距离跟圆的半径作比较. 但是椭圆不同了,而且旋转了 θ 角,让自然参数 t 也不直接对应于P点的方位了,只是稍微增加了一点点麻烦而已,仍然可以找出椭圆上的边界点,然后对比到椭圆圆心的距离. 这个方法太原始了,而且不易推广,显然不是我想说的.
看看math.stackexchange.com上人家怎么说吧.如何确定任一点是否在二次曲线内,我当然不是为了转一个链接而卖关子. 链接里的方法主要还是对方程为隐函数 g(x,y)=Ax2+Bxy+Cy2+Dx+Ey+F=0 形式的二次曲线,利用其二次型.但仅适用于二次曲线.
大部分情况下,参数方程更方便,而且两种形式之间互换本来也挺麻烦. 这里介绍基于椭圆的参数方程时,如何判断; 同时,把这个例子里的概念推广到更一般的高于二次的 C∞ 曲线，以及凸的和凹的多边形.
这个概念就是拓扑学中的winding number 或 卷绕数3.
很容易把任意曲线的参数方程:
{x=y=x(t)y(t)
改写成极坐标形式:
⎧⎩⎨⎪⎪⎪⎪⎪⎪ρ=θ=x(t)2+y(t)2−−−−−−−−−−√arctan(y(t)x(t)) 特殊的反正切,应该带卦限信息,并兼容x(t)=0
求导数之后(参考维基上的内容winding number 或 卷绕数):
dθ=x(t)y′(t)−y(t)x′(t)x(t)2+y(t)2dt
曲线 x(t)−y(t) 在某点 (x0,y0) 的卷绕数定义为:
w(x0,y0)=12π∫t1t0(x(t)−x0)y′(t)−(y(t)−y0)x′(t)(x(t)−x0)2+(y(t)−y0)2dt
对于前面提到的椭圆: t0=0,t1=2π ,**当 w(x0,y0) 为0时,点在椭圆外,为非0整数时在其内。
这个结论使得用数值积分计算结果作为判据非常合适,因为它对积分的数值计算有很大的允许误差; 从维基还是某本书上提到,winding number在量子力学里面又是一种”量子数”(稍后核实出处)。
用Mathematica代码做个试验:
curve={5Cos[t],3Sin[t]}; (*标准椭圆定义*)
{x[t],y[t]}=RotationTransform[Pi/]@curve; (*旋转30°的参数形式*)
innerQ[{x0_,y0_}]:=Round[NIntegrate[((x[t]-x0)D[y[t],t]-(y[t]-y0)D[x[t],t])/((x[t]-x0)^+(y[t]-y0)^),{t,0,2Pi}]//Pi]!= (*以卷绕数为判据的函数的定义*)
innerQ /@ {{3, 5}, {3, 2}, {0, 0}, {-4, 1}} (*测试几个点*)
{False, True, True, False}(*输出结果*)
太丑了! 换个表现方式!!
ClearAll["Global`*"]
curve={5Cos[t],3Sin[t]};
new=RotationTransform[Pi/]@curve;
x[t_]:=new[[]]
y[t_]:=new[[]]
innerQ[{x0_,y0_}]:=Round[NIntegrate[((x[t]-x0)D[y[t],t]-(y[t]-y0)D[x[t],t])/((x[t]-x0)^+(y[t]-y0)^),{t,0,2Pi},AccuracyGoal->]//Pi]!=
hc=ParametricPlot[u*RotationTransform[Pi/]@curve,{t,0,2Pi},{u,0,1},PlotRange->{{-5,5},{-4.5,4.5}},PlotPoints->,MeshFunctions->{#4&},Mesh->{{1.1}},MeshShading->{Yellow,Yellow}]/.Line[_]:>Sequence[];
SeedRandom[]
(points=RandomReal[{-4,4},{10,2}]);
txt=Text[Style[#,Red,],points[[#]],{1,0}]&/@Range[Length@points];
v=innerQ/@points;
Join[points,ConstantArray["->",{10,1}],Transpose@{Range[10]},Transpose@{v},]//MatrixForm
Show[hc,Graphics[{Blue,AbsolutePointSize[6],Point[#]&/@(1.15points)}],Epilog->{Red,txt},Frame->True,ImageSize->Large]
输出结果是这样滴:
| 编号 | 点的坐标 | 是否在曲线封闭区域之内 |
|---|---|---|
| 1 | (-3.99368, -3.47985) | 不在 |
| 2 | (3.91644, 3.75014) | 不在 |
| 3 | (-2.39307, 2.55617) | 不在 |
| 4 | (-3.28189, 3.76561) | 不在 |
| 5 | (-2.16072, 0.90002) | 在 |
| 6 | (-3.22547, 0.39084) | 在 |
| 7 | (-2.93962, -2.14134) | 在 |
| 8 | ( 2.20908, 0.40759 ) | 在 |
| 9 | (-3.53048, 3.68482) | 不在 |
| 10 | (-3.21401, -3.72518) | 不在 |
……………………………………………………………………………………………………………………………..
卷绕数的几何意义
这个部分希望能够增加一点直观生动的内容。
一个深入点的例子
前面的例子其实不能算简单，我开始后悔没有用”圆“了。
如果说圆和椭圆的情况，都还可以通过测量到”形心距离”之类的方法——即便有些小周折，这让基于winding number 或 卷绕数的方法还是有些不太自在的——”深入点的例子“必须给常规的潜在竞争关系的方法制造足够的麻烦，或者至少要让心目中想要突出的、中意的方法干得更加漂亮！
下面这个例子，一方面仍然保持图形以及参数方程本身在视觉上或表面看上去的足够简单，同时希望让那些初级的方法更加步履维艰，否则我写了这么多不是瞎忙活吗。这条曲线的参数方程是：
{x(t)=y(t)=9sin2t+5sin3t9cos2t−5cos3t
把它画出来是这样的，看上去是不是有些似曾相识？
在用凹多边形表述发现代码太长之后，我想尝试用曲线近似表达一颗五角星，结果就找到了它。它的边上的肌肤的线条都是 C∞ 光滑的，不像凹多边形那样有棱有角；但是形状有非常近似；所有的边界都可以用一个方程表达，看上去似乎有压缩代码的潜力。如果可以简洁高效地判断出内点外点，只须再作如下填充不就是五角星了吗：
代码悄悄地：
Grid@{{ParametricPlot[u {9 Sin[2 t]+5 Sin[3 t],9 Cos[2 t]-5 Cos[3 t]},{t,0,2 Pi},{u,0,1},MeshFunctions->{(*[email protected](#1^2+#2^2)*)#4&},Mesh->{{1}},PlotPoints->,MeshStyle->Red,MeshShading->{{Cyan,LightGreen}},ImageSize->Large]/.Line[x_]:>{Red,Line[x]},ParametricPlot[u {9 Sin[2 t]+5 Sin[3 t],9 Cos[2 t]-5 Cos[3 t]},{t,0,2 Pi},{u,0,1},MeshFunctions->{(*[email protected](#1^2+#2^2)*)#4&},Mesh->{{1}},PlotPoints->,MeshStyle->Red,MeshShading->{{Red,Red}},ImageSize->Large,Axes->False]/.Line[x_]:>{Red,Line[x]}}}
郁闷的是，对它进行判断的方法，跟椭圆的时候没有区别。所以，这部分就不讲了。连个偷着乐的表情都不能发，MarkDown太死板、不灵活。不过算法似乎可以趁机把MarkDown的流程图练习一下, 等着吧。
此外，对于这个例子为何比椭圆情况更进一步没有解释。
对winding number的解释,也缺乏维基那样生动的动态演示.
虽然不是必要的,作为练习还是可以慢慢加上来的.
向多边形的推广
惊艳的其中一部分在这里。 拓扑的确不一样。
卷绕数在五角星绘制中的应用实例
本来问题就是从绘制五星红旗开始的。为了压缩代码，希望找到一个新的判断内外点的算法。现在不论成功与否，贴代码吧。
我只把高级语言当作表达某些意思的工具，但是从来没有去研究关于这类工具的特殊的形式和技巧，无从着手，也不是特别感兴趣。所以，算法在此，进一步压缩只能翘首等别人了。
// NOTE: compile with g++ filename.cpp -std=c++11
#include <iostream>
#include <cmath>
#include <cstdlib>
#include <vector>
//#include <functional>
#define DIM1 600
#define DIM2 400
#define DM1 (DIM1-1)
#define DM2 (DIM2-1)
#define _sq(x) ((x)*(x)) // square
#define _cb(x) abs((x)*(x)*(x)) // absolute value of cube
#define _cr(x) (unsigned char)(pow((x),1.0/3.0)) // cube root
#define pi 3.14159
#define eps 2.02e-16
double arct(double x,double y){//可改写为宏放Rd函数中
return abs(x)>eps?(x>eps?atan(y/x):(y>?atan(y/x)+pi:atan(y/x)-pi)):(y>?pi/:-pi/);
}
bool inQ(int i,int j,int* px){//改为lambda表达式放GR/BL中
double art[],aro[];
for(int k=;k<;k++){
art[k]=arct((double)(px[k]-i),(double)(px[k+]-j));//
}
double total=;
for(int k=;k<;k++){
aro[k]=art[(k+)%]-art[k];
while(aro[k]<-pi) aro[k]+=pi;
while(aro[k]>pi) aro[k]-=pi;
total+=aro[k];
}
return abs(total)>;
};//即使修改了也只能作有限压缩，不能达到420字符水平
unsigned char GR(int,int);
unsigned char BL(int,int);
unsigned char RD(int i,int j){
// YOUR CODE HERE
return ;
}
unsigned char GR(int i,int j){
// YOUR CODE HERE
int p1[]={, , , , , , , , , },
p2[]={, , , , , , , , , },
p3[]={, , , , , , , , , },
p4[]={, , , , , , , , , },
p5[]={, , , , ,, , , , };
return inQ(i,j,p1)||inQ(i,j,p2)||inQ(i,j,p3)||inQ(i,j,p4)||inQ(i,j,p5)?:;
}
unsigned char BL(int i,int j){
// YOUR CODE HERE
return ;
}
void pixel_write(int,int);
FILE *fp;
int main(){
fp = fopen("MathPic.ppm","wb");
fprintf(fp, "P6\n%d %d\n255\n", DIM1, DIM2);
for(int j=;j<DIM2;j++)
for(int i=;i<DIM1;i++)
pixel_write(i,j);
fclose(fp);
system("pause");
return ;
}
void pixel_write(int i, int j){
static unsigned char color[];
color[] = RD(i,j)&;
color[] = GR(i,j)&;
color[] = BL(i,j)&;
fwrite(color, , , fp);
}
代码的输出结果ppm经过XnView或Convert转换成png格式之后是这样的：
foo bar foo
foo bar foo
Edit
【稍安勿躁，正在加班加点更新】
- “相濡以沫，不如相忘于江湖。” 我发现庄子就像读过鲁迅的《伤逝》一样，早在几千年前就建议好了该放下的时候要果断放下，可还是有人执迷不悟。——这就是古人的大智慧。“难道我们还不如古人吗？”(后面这个反问句在我高中一同学的作文里常用，每次想起他在班里当众读自己作文里这一句时气势汹汹的样子，我都想回答“是的。在很多方面。”这也是我瞧不起一些文科生和文科公知的原因: 立论完全罔顾事实,单纯依赖辩论技巧. 你他么这样自欺欺人能解决黎曼猜想还是哥德巴赫猜想?) ↩
- 我觉得刚才我写到这里没更新的时候,肯定有点进来的人看了很失望,因为实际上到此为止,除了画了个椭圆,我啥也没讲。会不会觉得博主不太靠谱？那一定是因为你自己更不靠谱！我觉得评价一个人靠谱与否的唯一标准是看此人是不是删过facebook ↩
- “拓扑学”听上去挺抽象，但我觉得数学之所以看上去让人望而却步，除了某些基础性的基本功的确需要积累和修炼之外，大部分因为老师的方法有问题。 在深入了解这个概念之前就知道了“判断是内点或外点”时这个概念该怎么用——也就是实际使用的方法是很简单的——，不过，我边写边完善的过程中发现有更多相关而有趣的不同表述，比我原先想象的精彩得多，需要自己细细琢磨，所以更新的速度可能会变慢。学习速度更快的不妨自己先看维基。 ↩