MATLAB绘制高等数学附录“地狱渐进难度”函数曲线

00.前言

       前几天有一个小伙伴问我怎么用MATLAB绘制一些函数的图像，我看了下，这不就是万恶的高等数学附录后面的图像内容吗？看到题目要求绘制的前六个图像，我的思绪又飘到了几年前被高等数学支配的恐惧中。

但是，这几个函数的确是锻炼MATLAB编码的好例子，因此特地写下这篇帖子，算是对那个u的反馈，对他信任我的MATLAB技术表示感谢。

先贴上此次任务的几个函数美貌让大家欣赏一下，见图1

图1：本次的目标函数群

01.绘制三次抛物线的图像

       分析：首先，我们需要注意，这是一个单调函数，而且在整个定义域内单调，这无疑是很简单的。

其次，我们需要考虑自变量  的取值范围。

不难分析出来  ，但是，必须注意，在计算机中很难将整个R域表达出来。

因此，我们只截取部分对称区间绘图。

在本题目中，取, 绘图的平滑度取决于自变量取值的步长大小。

因变量  随自变量  的变化而变化，这个不再赘述。

另外，请注意，这个函数中还有一个参数  ，虽然书中的图像并没有对其取值进行讨论，而且书中的  显然为正值。

但为了方便大家了解，本着严谨的态度，我们对其取值也进行讨论。

事实上，参数满足： ，但同样地，我们对其进行类似于自变量的简化处理，只考虑在某些特定条件下的取值。

      解：对于自变量  ，我们给步长 ，现对参数的取值进行讨论。

①当  时，我们只取  ，代码如下：

a = [-1.5 -2 -2.5 -5 -10];% 定义参数a的取值数组
x = -10:.5:10;% 定义自变量的取值范围，步长为0.5
for i =1:length(a)
    y = a(i) * x.^3; %定义函数关系式
    plot(x,y,'LineWidth',3)%绘制函数图像，定义线条宽度为3
    str{i} = ['a=',num2str(a(i))];%定义字符串数组str，用来存放随后需要用到的legend内容
    hold on;
end
legend(str)
title('当a < -1 时曲线的变化情况')
hold off

        代码执行结果如图2所示：

 图2：a<-1时的取值对函数的影响

        分析图2，可以发现： 的值越负，函数的开口越小，或者说函数越来越“陡峭”，但是函数的变化趋势没有发生变化。

细心的小伙伴可能已经注意到了，这与书中要求的函数图像不一致，这是不是你的代码有问题？别着急，我们接着往下分析。

②当 时，取  ，代码如下：

a = [-0.8 -0.5 -0.3 -0.1];
x = -10:.5:10;
for j = 1:length(a)
    y1 = a(j) * x.^3;
    plot(x,y1,'LineWidth',3)
    str_1{j} = ['a=',num2str(a(j))];
    hold on;
end
legend(str_1)
title('当 -1 < a <0时曲线的变化情况')
hold off

       代码运行结果如图3所示：

 图3：当 -1

        我们发现，此时的规律与①一致，两者其实可以合并，接下来讨论  取值为正的情况。

③当  时，取  进行讨论，话不多说，直接上代码：

a = [0.1 0.3 0.5 0.8];
x = -10:.5:10;
for j = 1:length(a)
    y2 = a(j) * x.^3;
    plot(x,y2,'LineWidth',3)
    str_2{j} = ['a=',num2str(a(j))];
    hold on;
end
legend(str_2)
title('当0 < a < 1时曲线的变化情况')
hold off

        代码运行结果如图4所示：

 图4：当 0

      此时，函数的整体变化趋势与之前相反，而且参数的值越大，函数的开口越小。

同样地，当 时，情况也是一样的，这里就不放了，感兴趣的小伙伴可以自己私下尝试。

      综上，如果想要实现题目要求的效果，只需要取参数为正就可以了，数值取多少什么的どうでもいい。

02.绘制半立方抛物线的图像

        分析：这个题目的自变量  取值范围已知，而且参数  的取值范围我们也可以直接给出： ，结合第一个函数分析过程，大家是不是觉得这样做很繁琐呢？其实这只是从逻辑的角度分析可能的状况，现在我们偷个懒，只考虑参数取整数的情况。

另外，如果从传统的逻辑角度分析这个函数，你将会发现一个问题，那就是：如果这是一个函数，那么必须满足：自变量取某一个值的时候，因变量都有唯一一个值与之一一对应。

但是，这个“函数“明显违背了这一个基本条件，怎么办呢？有两种解决思路：The first one，把这个关系式看成分段函数，然后分别对两半值域区间绘图即可。

The second one，将  看成自变量，将  看成因变量，这样，将满足函数定义。

现讨论第一种情况：

        The first one：很明显，以 0 为断点划分值域，那么两部分都是单调的，绘图将简便很多

①当  时，我们对关系式进行等价转换，可以得到 ，然后，上代码：

a = 1;% 为了简便起见，取参数a的值为1
x = 0:.5:10;%只定义自变量为正数

%这里相当于将因变量拆分成两个函数，然后在绘图时拼接即可。


y1 = -sqrt(a * x.^3);
y2 = sqrt(a * x.^3);

%绘图并且拼接
plot(x,y1,x,y2,'LineWidth',3)
hold on;
str_2 = ['a=',num2str(a)];
legend(str_2)
title('半立方抛物线y^2 = a*x^3图像')
hold off

        代码运行结果如图5(a)所示：

 图5(a)：半立方抛物线绘图结果

        第一种情况讨论完毕，我们来看看第二种情况。

       The second one，将表达式变形，得到  ，此时，我们不难发现，自变量的取值范围为：  ，另外的取值无任何变化。

话不多说，上代码：

a = 1;
y = -10:.5:10;
for i = 1:length(a) %这里都没必要做for循环
    x = (a*y.^2).^(1/3);
    plot(y,x)
    str = ['a = ',num2str(a)];
    hold on;
end
legend(str)
title('y^2 = x^3的图像')
hold off

        代码运行结果如图5(b)所示：

 图5(b)：半立方抛物线的运行结果

       不过这里有一个问题啊，跟第一种情况不一样，这里留一个小作业吧，怎么把这张图给它”摆正“？

03. 绘制概率曲线的图像

       分析：没什么好分析的，莽它就完了，直接上代码：

x = -3:.05:3;
y = exp(-x.^2);
plot(x,y)
title('y = e^-(x)^2')

        代码运行结果如图6所示：

 图6：概率函数的图像

高能预警！接下来的三个函数可不是跟你闹着玩儿的，搞不好的话，小则图画不出来，大则电脑卡死，玩家崩溃。

04.绘制箕舌线的图像

        分析：显然 ，自变量的取值范围为 ，为了获得和题目要求一致样貌的函数曲线，我们取。

需要注意的是，这个函数是由一个圆和上面的部分组成的。

可以考虑分部画图，需要指出的是，这个图可以直接用表达式作图，不需要利用极坐标或者参数方程完成。

好了，上代码：

a = 2;%取参数a的值为2，当然，你也可以取一个你喜欢的数值
x = -10:.1:10;%定义自变量的取值范围
y = 8*a^3./(x.^2+ 4*a^2);%定义函数公式

%画下面的那个圆
theta = 0:pi/20:2*pi;
xx = a*cos(theta);
yy = a+a*sin(theta);

%组装
plot(x,y,'b','LineWidth',3)
hold on
plot(xx,yy,'r--','LineWidth',1.5)
legend('箕舌线图像，a = 2')
axis equal

        代码运行结果如图7所示：

图7：箕舌线图像 

        配色稍微有点尴尬，如果你实在是看不下去可自行调整代码参数，绘制出属于你的amazing效果。

05.绘制蔓叶线的图像

        分析：这是一个隐式方程，可以直接用方程莽，但需要注意分母必须有意义的情况。

考虑画这个图的时候，让四五年没有碰高数课本的我着实死了好多脑细胞。

我们需要寻找这个方程的其它形式。

参考了EMBA乐园：蔓叶线漫谈https://m.sohu.com/a/505017388_121124316/?pvid=000115_3w_a        这一篇帖子。

我们发现，可以从常规的定义角度入手，推导这个曲线的极坐标方程。

现结合引文内容，与自己的思考，将推导过程再简述一遍。

虽然已经有前人做过这个工作了，但我们还是有必要研究一下这个过程。

放一张参考图图8作为说明：

 图8：分析过程专用图

       引用定义：

作一圆，取它的一条直径ON，不妨置于水平位置。

过N作ON的垂线l。

从点O引射线，交圆于点A，交垂线于点B。

以AB长度为半径，以点O为圆心作圆，与射线交于点M。

那么，随着射线的运动，点M就将描绘出一条曲线，这条曲线就是蔓叶线，

        为了方便起见，我们统一引用极坐标表达式表示其中的任何线条。

        对于直线  ，我们知道它的直角坐标表达式为： ，化为极坐标表达式为：。

如果不知道怎么换算，可以参考这篇文章：直线的极坐标方程https://zhuanlan.zhihu.com/p/55877974

     那么，在直线上的  点的极坐标就可以确定为：。

而  点在半径为a的圆上并且满足射线方程，那么其满足：，线段  的长度为：，其中，在中，我们可以很方便地利用直角三角形的角度和边的关系推导出点应该满足的关系式。

    即：。

根据引文中的转换，将其转换成直角坐标系的表达式，即：

。

话不多说，开始编写代码：

a = 2; %定义a的取值
x = 0:.05:4;%定义自变量的取值范围

ylim = get(gca,'Ylim');%获取当前窗口的Y轴数值


%定义主函数图像关系式
y1 = sqrt(x.^3./(2*a-x));
y2 = -sqrt(x.^3./(2*a-x));

%定义辅助圆参数
theta = -2*pi:pi/20:2*pi;
xx = a+a*cos(theta);
yy = a*sin(theta);

plot(x,y1,x,y2,'LIneWidth',3)
hold on
plot(xx,yy,[4,4],ylim,'LineWidth',1.5)
axis([0 10 -10 10])
grid on

    代码运行结果如图9所示:

图9： 蔓叶线绘图结果

 06.绘制笛卡尔叶形线的图像

       分析：这个堪称史诗级难度，不过幸好有参数方程，即：，有了参数方程就会好很多，其次，我们需要知道它的渐近线方程，即：。

上代码：

a = 1.5;
t = -5:.1:5;
x = 3*a.*t./(1+ t.^3);
y = 3*a.*t.^2./(1+ t.^3);
plot(x,y,'LineWidth',3)
hold on
grid on
x1 = -10:.001:10;
y1 = -a-x1;
plot(x1,y1,'--')
axis([-5 5 -5 5])

    运行结果如图10所示：

 图10：笛卡尔叶形线图像

         至于为什么中间那一段儿不封闭，我也不晓得咯~

07.后记

       从早上一直写到现在，花了一天时间了，终于写完了。

后面这几个着实废了一番功夫，但值得你花时间思考。

听说MATLAB还可以画出动态效果，这里就先不展示了，研究一段时间再说，如果大家感兴趣的话也可以自己尝试一下，很带劲儿。

 

---

       代码没有经过专门的优化，如果大家有更好的编写方式，欢迎留言讨论，我们下期再会，see you then~