使用Maxima进行极限微积分
我们在2011年2月刊中介绍了Maxima在代数运算和方程重排方面的应用。但这并非Maxima的全部技巧。本月，我将介绍Maxima如何帮助解决微分方程，但我将省略一些中间结果以节省空间。
许多科学领域都涉及弄清楚系统如何随时间变化以及导致这些变化的原因。当你开始研究变化，特别是变化率时，这本质上就是微积分。微积分和变化率也与图表上直线的斜率有关。当你绘制数据并找到描述它的方程时，可以通过对该方程求导数来找到直线的斜率。让我们看看一个下落的物体，看看理论对此有何说法。
你应该从研究如何求导数开始。假设你有方程
(%i1) f(x):= 2 + x^2;
2
(%o1) f(x) := 2 + x
你可以通过调用函数diff找到一阶导数，并提供要微分的方程以及要微分的变量。所以，你可以写成
(%i2) answer:diff(f(x),x); (%o2) 2x
Maxima也可以对表达式进行微分。如果你有几个方程，你可以用以下方法推导出它们的比率
(%i3) g(x):= x^(1/2);
(%i4) ratio_diff:diff(g(x)/f(x),x);
3/2
1 2 x
(%o4) ----------------- - -----------
2 2 2
2 sqrt(x) (x + 1) (x + 1)
这可能有点难以处理，所以你可能想要将其重构为更简洁的形式
(%i5) factor(ratio_diff);
2
3 x - 1
(%o5) - ------------------
2 2
2 sqrt(x) (x + 1)
Maxima也可以处理三角函数，但是有很多恒等式可以用来帮助简化包含三角函数的方程。默认情况下，除非你明确指出，使用特殊函数，否则Maxima不会尝试应用这些恒等式。例如，假设你有以下方程
(%i6) diff(sin(x)/(1 + cos(x)),x);
2
sin (x) cos(x)
(%o6) ------------- + ----------
2 cos(x) + 1
(cos(x) + 1)
(%i7) factor(%);
2 2
sin (x) + cos (x) + cos(x)
(%o7) --------------------------
2
(cos(x) + 1)
这仍然不是很简单。如果你然后应用函数trigsimp，你可以强制Maxima将三角简化规则应用于方程，并看看你会得到什么
(%i8) trigsimp(%);
1
(%o8) ----------
cos(x) + 1
你应该注意一些关于Maxima如何处理三角函数的重要的注意事项。首先是sin(x)^(-1)是正弦的倒数，而不是反正弦。要获得反正弦，你应该使用asin(x)。另一个是另一个三角简化函数，trigreduce。此函数用于通过使用多倍角公式来降低三角函数的幂。例如
(%i9) trigsimp(cos(x)^2 + 2*sin(x)^2);
2
(%o9) sin (x) + 1
(%i10) trigreduce(cos(x)^2 + 2*sin(x)^2);
cos(2 x) + 1 1 cos(2 x)
(%o10) ------------ + 2 (- - --------)
2 2 2
这看起来可能没有trigsimp的结果简单，但它是方程的一种更简单的形式，可以与其他函数一起使用，例如积分。
当求导数时，Maxima可以应用链式法则。假设你有方程
(%i11) f(x):= x^3);
3
(%o11) f(x) := x
(%i12) depends(x,u)$
(%i13) diff(f(x),u);
2 dx
(%o13) 3 x --
du
%i12行的代码使用了一个新函数depends。这是一种告诉Maxima x是u的函数的方法，而无需显式定义描述这种关系的函数。如果你稍后决定要为此关系定义一个实际的方程，你可以使用
(%i14) remove([x,u],dependency);
(%o14) done
(%i15) x:sin(u);
(%o15) sin(u)
(%i16) diff(f(x),u);
2
(%o16) 3 cos(u) sin (u)
同样地，Maxima可以处理隐式微分。假设你有方程 x^2 + y^2 = 25，并且你想求 dy/dx。你需要使用我刚才提到的depends函数来处理这个问题
(%i17) eqn := x^2 + y^2 = 25;
2 2
(%o17) y + x = 25
(%i18) depends(y,x);
(%o18) [y(x)]
(%i19) deriv_of_eqn:diff(eqn,x);
dy
(%o19) 2 y -- + 2 x = 0
dx
(%i20) solve(deriv_of_eqn,'diff(y,x));
dy x
(%o20) [-- = - -]
dx y
微积分的另一方面是积分。在Maxima中执行积分的基本函数称为integrate。此函数可以执行定积分和不定积分。不定积分是你在学校可能学到的积分的符号形式。例如
(%i21) integrate(x^2,x);
3
x
(%o21) -
3
定积分实际上是在一个区间上求值的。这种形式的积分可以可视化为由你正在积分的方程定义的曲线下的面积。要进行定积分，只需添加两个参数，给出区间的起点和终点
(%i22) integrate(x^2,x,0,1);
1
(%o22) -
3
将所有这些技术结合在一起，你可以求解关于给定变量的微分方程——例如，求解关于y的dy/dx = f(x)。你可以通过执行所有必需的代数和微积分来做到这一点，但你真的不需要。Maxima有一个非常强大的函数ode2，它可以一步完成。从你的方程开始
(%i23) eq: 'diff(y,x) = sqrt(1/x^2 - 1/x^3);
dy 1 1
(%o23) -- = sqrt(-- - --)
dx 2 3
x x
(%i24) ode2(eq,y,x);
2
2 2 sqrt(x - x)
(%o24) y = log(2 sqrt(x - x) + 2 x - 1) - ------------- + %c
x
这个函数调用执行积分和求解步骤，并给你微分方程的最终答案。
假设你正在做一个实验，扔下一枚硬币并计时它下落的时间。你如何知道你测量的时间实际上是否有意义？让我们从最基本的定律开始：力 = 质量 * 加速度。
硬币的质量是常数，所以暂时忽略它。力是由于重力，将硬币拉向地面，加速度描述了硬币由于这种力而产生的运动。重力是常数，至少在地球上是这样，并且它与质量呈线性关系，所以你可以将力定义为
(%i1) force: mass * g; (%o1) g mass
加速度也是一个常数，因为质量和力都是常数。加速度只是速度的变化率，速度是位置的变化率，所以你可以这样设置
(%i2) depends(y,t);
(%o2) [y(t)]
(%i3) acceleration: 'diff('diff(y,t),t);
2
d y
(%o3) ---
2
dt
将它们放在一起，你得到
(%i4) eq_of_motion: force = mass * acceleration;
2
d y
(%o4) g mass = mass ---
2
dt
(%i5) solve(eq_of_motion, y);
2
d y
(%o5) [--- = g]
2
dt
你可以立即看到物体下落的速度根本不取决于质量。伽利略是对的！下一步是做一些积分，看看你最终得到什么
(%i6) integrate(%,t);
dy
(%o6) [-- = g t + %c1]
dt
在这一步，你将能够找出在时间t的速度（dy/dt）。附加项%c1是积分常数。在这种情况下，你可以看到它代表你的便士的初始速度。再进行一轮积分得到
(%i7) integrate(%,t);
/ 2
[ dy g t
(%o7) [I -- dt = ---- + %c1 t + %c2]
] dt 2
/
现在你可以找到你的硬币在任何时间t的位置y。再次，引入一个新的积分常数%c2。在这种情况下，你可以看到这代表你的硬币的起始高度。但这不是你正在测量的。你正在测量硬币下落给定距离所需的时间。所以你需要做一些重新排列。因为你正在扔下你的硬币，你知道起始速度是0（即%c1=0）。你可以稍微重写一下，使其更清晰
(%i8) eqn: y = (g * t^2)/2 + y0;
2
g t
(%o8) y = y0 + ----
2
(%i9) solve(eqn,t);
y y0 y y0
(%o9) [t = - sqrt(2) sqrt(- - --), t = sqrt(2) sqrt(- - --)]
g g g g
就这样。现在你有一个关于时间的方程，给定你的硬币下落的高度。有了这个理论关系，你可以检查重力是否在你的本地实验室中正常工作。如果不是，你应该立即联系诺贝尔委员会。
这仅仅触及了Maxima在处理微积分和微分方程方面的能力的表面，但希望本文能为你提供一个起点。祝你积分愉快。
