Matlab绘制双纽线、莫比乌斯环，双纽线公式、重积分

1. 引言：

今天在做一道二重积分时遇到了这样一个问题，看题：

设 $D=\{(x,y)|(x^2+y^2{)}^2\le 4(x^2-y^2)\}$ ，则区域D的面积为?

当时我还不知道这是双纽线的公式，就直接这样做了：

算出来面积怎么等于0呢？

原因就在于没有准确考虑$r和\theta$的取值；如果知道这个积分区域就是双纽线的内部，会更容易理解计算过程。

2. 双纽线公式：

伯努利双纽线：

$r^2=a^{2}sin2\theta$ 或 $r^2=a^{2}cos2\theta$

Matlab绘制代码：

clc,clear,close all;
theta=0:0.001:2*pi;
subplot(121)
r=sqrt(sin(2*theta));  
polarplot(theta,r,'r','linewidth',1.5);
title('r^2=sin2\theta')
subplot(122)
r1=sqrt(cos(2*theta));   
polarplot(theta,r1,'b--','linewidth',1.5);
title('r^2=cos2\theta')
1
2
3
4
5
6
7
8
9
10

将公式稍微修改一下，还可以得到其他有趣的图形：

clc,clear,close all;
theta=0:0.001:2*pi;
subplot(221)
r=sin(3*theta);  
polarplot(theta,r,'r','linewidth',1.5);
title('r=sin3\theta')
subplot(222)
r1=cos(3*theta);   
polarplot(theta,r1,'linewidth',1.5);
title('r=cos3\theta')
subplot(223)
r=sin(2*theta);  
polarplot(theta,r,'g','linewidth',1.5);
title('r=sin2\theta')
subplot(224)
r1=cos(2*theta);   
polarplot(theta,r1,'c','linewidth',1.5);
title('r=cos2\theta')

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

这里看子图中最经典的伯努利双纽线，是不是有点似曾相识的感觉呢？

诶对了，其实伯努利双纽线非常像一个二维的对称式双螺旋-莫比乌斯环，类似于极限的符号，本身也是永恒的象征；

在很多软件的logo、影视造型里也用到了莫比乌斯环，和双纽线的效果，比如…童年熟悉的梦比优斯奥特曼（实际上名字就是源于莫比乌斯带）

以及前段时间很火的美剧《洛基》中，时间管理局的侦探角色-莫比乌斯：

当然还有复联4中唐尼随手一个建模成功就摆出来了莫比乌斯环，作为穿越量子宇宙回到过去时间点…（此处省略三千个…）

甚至Zhi乎上还有各路神仙对此展开了详细的原理解释，有空闲时间的同学可以去围观：

3. 双纽线二重积分计算：

回过头来看看引言中的那道题，

当 $\theta$ 转到$\pi /4$时半径 $r$ 已经为0，故外层积分上限不应直接对 $\theta$ 取$(0,\pi /2)$或$(0,2\pi )$，那样是错误的，且根据对称性，可以通过以下方法计算出D的面积：

再看一道例题：

设区域D由双纽线 $(x^2+y^2{)}^2=x^2-y^2$ 围成，则 $I=\int {\int }_D\frac{x^2+2y}{x^2+y^2}dxdy=$ ？

解：

首先根据对称性，被积函数分子中的2y是关于积分区域 $x$ 轴上半部分和 $x$ 轴下半部分奇对称的，故带 $2y$ 项积分为0，直接忽略，剩下的转极坐标，并化简后利用华里士公式计算：

4. Matlab绘制莫比乌斯环：

莫比乌斯带由德国数学家莫比乌斯（Mobius，1790～1868）和约翰·李斯丁于1858年发现。就是把一根纸条扭转180°后，两头再粘接起来做成的纸带圈，具有魔术般的性质。
普通纸带具有两个面（即双侧曲面），一个正面，一个反面，两个面可以涂成不同的颜色；而这样的纸带只有一个面（即单侧曲面），一只小虫可以爬遍整个曲面而不必跨过它的边缘。这种纸带被称为“莫比乌斯带”（也就是说，它的曲面从两个减少到只有一个）。

Matlab代码：

clc,clear;
k = 10000;
u = 0:pi/k:2*pi;
v = [-1 -0.5 0 0.5 1];
for j = 1:length(v)
for i = 1:length(u)
x(i,j) = (1+v(j)*cos(u(i)/2)/2)*cos(u(i));
y(i,j) = (1+v(j)*cos(u(i)/2)/2)*sin(u(i));
z(i,j) = v(j)*sin(u(i)/2)/2;
end
end
axis vis3d  %3维坐标系
mesh(x,y,z)
title("单螺旋莫比乌斯环")
c = autumn(256);
colormap(c);
shading interp;
for i=1:36
camorbit(10,0,'data',[0,0,1])%[0 0 1]表示按z轴旋转。36*10=360表示旋转一周
M=getframe(gcf);
nn=frame2im(M);
[nn,cm]=rgb2ind(nn,256);
if i==1
imwrite(nn,cm,'out.gif','gif','LoopCount',inf,'DelayTime',0.1);%说明loopcount只是在i==1的时候才有用
else
imwrite(nn,cm,'out.gif','gif','WriteMode','append','DelayTime',0.1)%当i>=2的时候loopcount不起作用
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

感兴趣的同学可以尝试用Matlab绘制出双螺旋的莫比乌斯环动画！