常微分方程组的SAGE快速入门

这 Sage 快速入门教程是为MAA预备研讨会“SAGE：在本科生中使用开源数学软件”开发的(由美国国家科学基金会提供，截止日期0817071)。它是在知识共享署名-Share Alike 3.0许可证下授权的 (CC BY-SA )。

求解微分方程式是精确方法和数值方法的结合，因此是用计算机探索的好地方。我们已经在微积分教程中看到了一个这样的例子，值得回顾一下。

基本的符号技巧

sage: y = function('y')(x)
sage: de = diff(y,x) + y -2
sage: h = desolve(de, y)

暂时忘掉绘图，请注意，要象征性地解微分方程，需要做三件事：

  • 抽象函数 y ；

  • 一个微分方程式，我们把它放在单独的一行里；

  • 实际的 d n不同 e 方程式 solve 命令(首字母缩写为粗体 desolve )。

sage: show(expand(h))
\[C e^{\Left(-x\Right)}+2\]

因为我们没有指定任何初始条件，所以Sage(来自Maxima)放入一个参数。如果我们想要设置一个初始条件，我们使用 ics (适用于 i 初始 c 条件s )。例如，我们设置 ics=[0,3] 要指定何时 x=0 ， y=3 。

sage: h = desolve(de, y, ics=[0,3]); h
(2*e^x + 1)*e^(-x)

当然，我们已经注意到，我们可以用斜率场来描绘这一切。

sage: y = var('y') # Needed so we can plot
sage: Plot1=plot_slope_field(2-y,(x,0,3),(y,0,5))
sage: Plot2=plot(h,x,0,3)
sage: Plot1+Plot2
Graphics object consisting of 2 graphics primitives

备注

关于符号函数与符号变量：

  • 如果你想让 y 抽象函数，而不是变量，你必须单独做这件事。一个微分方程式需要一个 function 但谋划需要一个 var 。

  • 另一种选择是让 z 是垂直轴变量的名称。

  • 无论哪种方式，都会有一些东西不得不做出让步，因为通常来讲，我们对待这些事情 y 既是一个变量又是一个函数，这要用计算机来完成要复杂得多。

在Sage中还有许多其他的微分方程式工具。我们不能在快速入门中涵盖这一点的所有变体，但文档对符号解算器很有帮助。

sage: desolvers?

例如，Maxima可以做系统，也可以使用拉普拉斯变换，为了便于使用，我们包含了这些变换的版本。

在所有的微分方程解例程中，注意语法是很重要的！在下面的示例中，我们将微分方程放在命令的主体中，并且必须指定 f 是 d 从属 var 可用 (dvar )，并给出初始条件 \(f(0)=1\) 和 \(f'(0)=2\) ，它给出了示例中的最后一个列表。

sage: f=function('f')(x)
sage: desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2])
x*e^x + e^x

sage: g(x)=x*e^x+e^x
sage: derivative(g,x,2)-2*derivative(g,x)+g
x |--> 0

数值方法和幂函数级数方法

也有数值方法。

例如，上面的选项之一是 desolve_rk4 。这是一种四阶龙格-库塔方法，并返回适当的(数值)输出。在这里，我们 must 给出因变量 and 初始条件。

sage: y = function('y')(x)
sage: de = diff(y,x) + y -2
sage: h = desolve_rk4(de, y, step=.05, ics=[0,3])

将其与原始的象征性解决方案进行比较可能会很有趣。我们使用 points 来自高级打印教程的命令。

sage: h1 = desolve(de, y, ics=[0,3])
sage: plot(h1,(x,0,5),color='red')+points(h)
Graphics object consisting of 2 graphics primitives

这里的数字例程的主要用途是教学性质的。

对于更高级的数值程序，我们主要使用GNU科学库。使用这种方法稍微复杂一些，但提供了丰富的选择。

sage: ode_solver?

我们甚至可以做POWER系列解决方案。为了做到这一点，我们必须首先定义一个特殊的 power series ring 包括精确度。

sage: R.<t> = PowerSeriesRing(QQ, default_prec=10)
sage: a = -1 + 0*t
sage: b = 2 + 0*t
sage: h=a.solve_linear_de(b=b,f0=3,prec=10)

这个POWER系列解决方案在一段时间内相当不错！

sage: h = h.polynomial()
sage: plot(h,-2,5)+plot(2+e^-x,(x,-2,5),color='red',linestyle=':',thickness=3)
Graphics object consisting of 2 graphics primitives

这只是一个介绍；在其他地方有很多关于使用Sage的微分方程的资源，包括David Joyner的一本书，他为Sage编写了大部分原始代码包装Maxima来实现这一点。