matlab在解析几何上的简单应用？

写点废话

　　这学期我教本科生《解析几何》课程，我的研究方向和几何差别很大，但是教课过程中觉得这门课还是蛮有意思的，这已经是我第二次教这门课了，这一次明显感觉比上一次轻松。

　　说到解析几何这门课程，我就不得不提一下我自己本科的时候我们班这门课的“趣（悲）事（剧）”，我们当时是两个行政班一起上这门课，两个班加一起59人，这门课挂了29人，补考只过了一位同学。当时我们学校的论坛这件事直接上了“今日十大”版块。这位解析几何老师同时教二年级的《数值分析》课程，听说那门课挂的更惨，一时间这位老师在论坛上被骂的很惨。等到大二的时候，我们的《数值分析》还是这位老师教，大家都觉得完了，但是我们的《数值分析》课程几乎没有人挂科，教完我们之后，这位老师就转岗到行政部门了，当然，他的转岗也许是因为学校的要求（没有博士学位的教师不能上讲台授课），也许是因为对我们失望。我和这位老师私下没有过多的交集，他是1975年的，比我大了十多岁，他讲课思路很清晰，我很喜欢听他的课，我记得他开了门公选课叫《数学的源与流》，我去蹭过好多次。他上课还是比较严肃的，他我们的时候已经34岁了，我离那个年纪还有几年，这也许是我和他风格不一样的原因，也许过几年我上课也会比现在严肃一些。听说他现在在学校教务处已经是处级干部了，希望以后有机会回去还能向他讨教。

　　这次教这门课，我尝试给学生推荐matlab软件，这个软件确实功能强大，用微信公众号上的话说，“除了生孩子，matlab啥都能干”。周一的课内容比较少，所以我用了一节课多一点的时间给大家简单介绍了这个软件，主要是告诉大家matlab在解析几何上的一些作用。下面我按照我课上讲解的顺序大概整理一下。今年这门课我教三个大班，可能在每个班讲的内容略微有点区别，但是总体差不多。

高等代数和数学分析上的简单应用

　　首先，我给大家介绍了一点matlab在另外两门基础课《高等代数》和《数学分析》课程上的应用，考虑到大家目前只学了刚开始的内容，所以我只介绍了相应的内容。

高等代数——多项式计算

　　我在给大家讲解析几何课程《1.5 标架与坐标》的时候就给大家提过，计算机怎么保存一个多项式方便？通过讨论，我们大概知道了，计算机保存多项式只需要保存系数就可以，它们的系数就可以看成它们的坐标。这里matlab对多项式的处理都是基于它的系数进行，通过下面的命令和运行结果相信大家能明白一个系数怎么和一个多项式对应。

p1=[1 2 3 0 4];
poly2str(p1,\'x\')%将系数转换成多项式，字母用\'x\'

　　它的运行结果是

x^4 + 2 x^3 + 3 x^2 + 4

　　两个多项式的加减自然也就可以利用两个向量的加减完成，这里需要注意的是，如果次数不一样，得把低次的那个向量高位补0。加减无所谓，这两个运算也是比较简单的，那么多项式乘法该如何进行呢？下面就是体现matlab强大功能的一个开始。话不多说，上代码。

p1=[1 2 1];
p2=[1 1];
p3=conv(p1,p2);
poly2str(p1,\'x\')
poly2str(p2,\'x\')
poly2str(p3,\'x\')

　　运行结果是

x^2 + 2 x + 1
x + 1
x^3 + 3 x^2 + 3 x + 1

　　大家可以看到，多项式乘法，一个conv函数就可以搞定。可能高等代数最开始的时候乘法不怎么让大家难过，除法才是让人头疼的，matlab怎么做多项式的除法呢？我们知道，除法除完有商和余式，其代码如下：

p1=[1 3 3 1];
p2=[1 2];
[a r]=deconv(p1,p2);
poly2str(a,\'x\')
poly2str(r,\'x\')

　　运行结果如下：

  x^2 +  x + 1
  -1

　　这两个结果分别对应上面的a和r，大家可以自行验算一下。事实上，matlab还可以直接计算复数范围内多项式的根（零点），而且命令简单的令人发指：

p=[1 3 3 1];
roots(p)

　　运行结果是：

  -1.0000 + 0.0000i
  -1.0000 - 0.0000i
  -1.0000 + 0.0000i

　　我们可以看到，输出了3重根。关于多项式的相关计算我们就介绍这么多。

高等代数——行列式计算

　　虽然大家高等代数才学行列式不久，但是我在解析几何的课程上已经给大家讲过行列式了，这是附录部分的内容，也是解析几何课程的工具，大家应该能感觉到计算起来很麻烦，但是对于matlab也就是一个命令的事。

A=[1 2 3;
   2 3 4;
   1 1 7];
det(A)

　　运行结果

-6

　　这只是小试牛刀，matlab其实就是“Matrix+Laborary”的缩写，它处理矩阵的工具很多，什么特征值、特征向量、逆矩阵都是一个命令的事。

数学分析——求极限

　　下面我们介绍一点matlab在数学分析上的应用，还是考虑到大家目前只学了求极限，所以我们先介绍求极限的相关matlab命令。这里顺便提一下，自然对数底数e在matlab中是不认识的，但是matlab中对圆周率$\pi$是固化的，就是pi，命令和结果截图如下：

　　在matlab中，如果遇到要用e表示的数，它会用exp($\cdot$)显示，比如说咱们学的重要极限：

syms x
limit(\'(1+1/x)^x\',x,inf)

　　运行结果：

exp(1)

　　它比我们想象的要聪明，例如极限$\lim\limits_{x\rightarrow\infty}(1+\dfrac{1}{ax})^{x}$，我觉得并不是所有大一学生都能轻松算出它的结果是$e^{1/a}$，但是matlab可以：

syms a x
limit(\'(1+1/(a*x))^x\',x,inf)

　　运行结果：

exp(1/a)

　　我不想解释上面那个求极限的命令，我再给大家写一个求$\lim\limits_{x\rightarrow2}\dfrac{x^2-4}{x-2}$的命令：

syms x
limit(\'(x^2-4)/(x-2)\',x,2)

　　运行结果：

4

　　多的我也不说了，大家可以自行探索。

数学分析——函数作图

　　对于上面那些极限大家可以求出来，但是我们有时候会遇到一些很难想象的极限，比如$\lim\limits_{x\rightarrow0}x\sin\frac{1}{x}=0$，这是无穷小乘以有界变量，自然是无穷小，那么它是怎么收敛的呢？再比如，$\lim\limits_{x\rightarrow0}\sin\frac{1}{x}$的极限不存在，那么在$x\rightarrow0$时，这个函数图像又是怎么样的呢？这里，我们可以利用matlab做函数图像帮我们理解上面的问题。

　　不知道大家有没有想过，计算机怎么作图的？其实，计算机作图的本质是一个个的把点打到屏幕上，matlab也是这样，matlab里面绘制平面曲线的命令是plot函数，我们可以看个例子：

plot(2,3,\'b*\');%b表示蓝色，*表示线条形状

　　结果如下：

　　下面我们考虑$y=\sin(x)$的五点法作图，命令如下：

x=0:0.5*pi:2*pi;
y=sin(x);
plot(x,y,\'b*\');

　　结果如下：

　　我们只描出了五个点，中学老师告诉我们，用平滑的曲线把这几个点连起来就是正弦函数的图像了，但是matlab这时候就不那么聪明了，它不知道什么叫做平滑，它是直性子，你让它把点连起来，它会用直线：

x=0:0.5*pi:2*pi;
y=sin(x);
plot(x,y,\'b\');%默认线条是直线,即用直线连接相邻两点

　　结果如下：

　　很明显，这不是正弦函数的图像，我们尝试多弄几个点，也就是x的步长取小一些：

x=0:0.5:2*pi;
y=sin(x);
plot(x,y,\'b*\');

　　结果如下：

　　我们再用直线连起来：

　　大家可以想想，教材上那些比较光滑的曲线都是怎么画出来的，其实就是讲步长取得足够小，在这种情况下曲线和直线是一样的，我们可以将步长取为0.1，得到的效果如下：

　　看上去这个图很光滑了，但是实际上，计算机还是一个个点打上去，然后用直线连接相邻的点，其点状图如下：

　　计算机就这样“骗”了我们，给我们呈现了一条“光滑的曲线”。大家可以尝试绘制刚才提及的两个极限，这里我就不给出代码了。$y=x\sin\frac{1}{x}$，$x\in[-0.1,0.1]$的图像如下：

$y=\sin\frac{1}{x}$，$x\in[-0.1,0.1]$的图像如下：

　　我们可以看到，在$x\rightarrow0$时，前者确实趋于0，后者一直振荡，极限不存在。

解析几何上的简单应用

　　做了那么多铺垫，下面可以讲讲它在解析几何这门课程上的应用了，终于等到你，还好我没放弃。

点积和叉乘

　　第一章相对麻烦一点的计算就是点积和叉乘，matlab可以很轻松的完成基于坐标（右手空间直角坐标系）的点积和叉乘，话不多少，上代码和结果：

　　我想我不需要过多的解释了，大家可以想想，混合积该怎么处理，实际上也就是dot(cross(a,b),c)即可。

平面曲线绘制

　　前面我们讲了如何绘制函数图像，这一部分显得有点多余，但事实并不是这样。（一元）函数图像可以看成平面曲线，但是平面曲线不一定可以看成某个函数的图像，因此完全靠前面的方法是不行的。比如说，我们如何绘制单位圆$x^2+y^2=1$在平面上的曲线呢？可能会有人想到，将它转化成$y=\sqrt{1-x^2}$，这个当然可以，不过你只能绘制一个半圆，还有一半你得再写一个$y=-\sqrt{1-x^2}$再绘制函数图像，两部分拼起来可以得到一个完整的圆，其代码和结果如下，为了区分，我用不同的颜色表示上下两部分：

x=-1:0.01:1;
y1=sqrt(1-x.^2);
y2=-sqrt(1-x.^2);
plot(x,y1,\'b\',x,y2,\'r\');
axis equal%这一句是为了让整张图是正方形

　　很多曲线的普通方程是很难得到的，比如咱们学的摆线，就算勉强得到了普通方程，将它写成$y=f(x)$的形式也是很困难的，因此，我们绘制曲线的时候，常常借助参数方程，比如说椭圆的曲线，我们可以利用其参数方程：

a=5;b=3;
theta=0:0.01:2*pi;
x=a*cos(theta);
y=b*sin(theta);
plot(x,y,\'b\');
axis equal

　　我估计这个代码我不需要解释太多，你们写出椭圆的方程对比代码就能看明白。下面我们可以看看我们教材上的几个例子，我给出摆线的代码，并解释一下，其余的线条比如内摆线、外摆线、渐伸线等我就给出结果，大家可以模仿着自己谢谢代码。

clear all;
clc;
a=10;
theta=0:0.02:2*pi;
axis equal;
axis([-a a*max(theta)+a 0 30]);
F=getframe(gcf);
m=frame2im(F);
[m,cm]=rgb2ind(m,256);
imwrite(m,cm,\'baixian.gif\',\'gif\');
for k=1:length(theta)
x=a*(theta(1:k)-sin(theta(1:k)));
y=a*(1-cos(theta(1:k)));
theta0=0:0.01:2*pi;
x1=a*theta(k)+a*cos(theta0);
y1=a+a*sin(theta0);
plot(x1,y1,\'b\',x,y,\'r.\');
axis equal;
axis([-a a*max(theta)+a 0 30]);
%%%%生成gif格式图片-start
F=getframe(gcf);
m=frame2im(F);
[m,cm]=rgb2ind(m,256);
imwrite(m,cm,\'baixian.gif\',\'gif\',\'WriteMode\',\'append\',\'DelayTime\',0.00001);
%%%%生成gif格式图片
end

　　上面就是程序运行的结果和代码，我简单解释一下，这里的动画每一帧都有两条曲线，一条是我们的摆线，另一条是滚动的圆。滚动的圆实际上就是圆心不断变化的圆，在我们的程序里是用x1和y1来表示的，大家可以看到，x1和y1就是圆心在不断变化，半径始终为$a$的圆的参数方程，另外，咱们的摆线的参数方程教材上有，大家可以自己翻翻书，再对比咱们的代码。最后有几行代码是保存成gif格式的图片的，这个是通用的。下面我们可以看看内摆线，按照教材上写的，如果$a=4b$，结果应该如下：

　　如果我们设置$a=10b$，效果如下：

　　外摆线是课后习题，我们给出$a=b$的“心脏线”的效果如下：

　　这两个，我希望大家可以自己动手尝试，模仿我上面给出的摆线的程序，自己谢谢内外摆线的程序，区别就是，内外摆线的程序每一帧都是要绘制三个曲线，一个定圆，一个动圆，一个就是我们的内外摆线，再根据教材提供的相关参数方程即可得到此代码。希望大家明白，计算机可以帮我们计算，但是不能帮我们建模！

空间曲面绘制

　　上面是平面上的曲线，内容相对简单些，下面我们考虑空间曲面的绘制，空间曲面的绘制相对而言多了一步将数据利用meshgrid函数网格化，一般来说，空间曲面可以看成是二元函数的图像，或者看成有两个参数的参数方程。我个人比较喜欢用参数方程绘制空间曲面，例如球面，根据其参数方程，代码和效果如下：

a=5;%球的半径
phi=0:0.01:2*pi;
theta=-0.5*pi:0.01:0.5*pi;
[phi theta]=meshgrid(phi,theta);
x=a*cos(theta).*cos(phi);
y=a*cos(theta).*sin(phi);
z=a*sin(theta);
mesh(x,y,z);
axis equal

　　类似的，柱面的绘制根据其参数方程代码和效果如下：

a=5;
phi=0:0.01:2*pi;
u=-10:0.1:80;
[phi u]=meshgrid(phi,u);
x=a*cos(phi);
y=a*sin(phi);
z=u;
mesh(x,y,z);

空间曲线绘制

　　空间曲线的绘制和平面曲线的绘制比较类似，都是打点，用的函数是plot3函数，如果大家对前面的plot函数能理解的话，这里代码就很容易理解了，我们可以考虑教材上的例子圆柱螺旋线，根据其参数方程，代码如下：

a=5;
t=0:0.01:12;
omega=3;
b=2;
x=a*cos(omega*t);
y=a*sin(omega*t);
z=b*omega*t;
plot3(x,y,z,\'r.\');

杂例——螺旋线动画

　　下面我们将螺旋线和圆柱面在一张图上显示，并且让螺旋线上的点逐个打出，生成gif图片，matlab中hold on命令可以保持原来的图像不被新的图像覆盖，代码如下，大家可自己运行观察效果：

a=5;
phi=0:0.01:2*pi;
u=0:0.01:8;
[phi u]=meshgrid(phi,u);
x=a*cos(phi);
y=a*sin(phi);
z=u;
mesh(x,y,z);
hold on;
t=0:0.01:12;
omega=3;
b=0.2;
x=a*cos(omega*t);
y=a*sin(omega*t);
z=b*omega*t;
for k=1:length(x)
    plot3(x(k),y(k),z(k),\'r.\');
    hold on;
    pause(0.00001)
end

杂例——心形线

　　数学是很浪漫的，讲了这么多曲线，就出现了一个心脏线，看着还很别扭，实际上，很多方程都能产生心形曲线，这里给出两个简单的例子，我们可以利用ezplot函数绘制隐函数确定的曲线，直接得到结果：

1. $x^2+(y-(x^2)^{\frac{1}{3}})^2=9$

ezplot(\'x^2+(y-(x^2)^(1/3))^2=9\');

1. $-x^2y^3+(x^2+y^2-1)^3=0$

ezplot(\'-x^2*y^3+(x^2+y^2-1)^3=0\',[-1.5,1.5]);

1. $17x^2-16|x|y+17y^2=200$

ezplot(\'17*x.^2-16*abs(x).*y+17*y.^2=200\');

杂例——心形曲面

　　平面上的心大家可能觉得不过瘾，我们再提供一个心形曲面的方程吧：

$(x^2+ \dfrac{9}{4}y^2 + z^2 - 1)^3 - x^2z^3 - \dfrac{9}{80}y^2z^3=0$

　　其效果如下：

　　如果你觉得不好看，可以适当调整，结果如下：

　　这两个心形曲面的代码依次为：

f=@(x,y,z)(x.^2+ (9./4).*y.^2 + z.^2 - 1).^3 - x.^2.*z.^3 - (9./80).*y.^2.*z.^3;
[x,y,z]=meshgrid(linspace(-1.5,1.5));
val=f(x,y,z);
isosurface(x,y,z,val,0); 
axis equal;view(3);colormap([1 0.2 0.2]);

　　和

f=@(x,y,z)(x.^2+ (9./4).*y.^2 + z.^2 - 1).^3 - x.^2.*z.^3 - (9./80).*y.^2.*z.^3;
[x,y,z]=meshgrid(linspace(-3,3));
val=f(x,y,z);
[p,v]=isosurface(x,y,z,val,0);
patch(\'faces\',p,\'vertices\',v,\'facevertexcdata\',jet(size(v,1)),\'facecolor\',\'w\',\'edgecolor\',\'flat\');
view(3);grid on;axis equal;

2019.11.12-2019.11.13 王飞