基本的代数和微积分

Sage可以执行许多基本的代数和微积分运算：如方程求解,微分,积分和Laplace变换。更多例子参见Sage的构成。

注意在下面的例子中,函数中出现的变量需要事先使用var(...)来定义。比如：

  1. sage: u = var('u')
  2. sage: diff(sin(u), u)
  3. cos(u)

如果出现了NameError错误,检查你是否有拼写错误,或者是忘记了用var(...)来定义这个变量

解方程

精确求解方程

solve函数用于解方程。要使用它,先要指定变量,然后将方程(或方程组)以及要求解的变量作为参数传给solve. 

  1. sage: x = var('x')
  2. sage: solve(x^2 + 3*x + 2, x)
  3. [x == -2, x == -1]

可以求解单变量的方程：

  1. sage: x, b, c = var('x b c')
  2. sage: solve([x^2 + b*x + c == 0],x)
  3. [x == -1/2*b - 1/2*sqrt(b^2 - 4*c), x == -1/2*b + 1/2*sqrt(b^2 - 4*c)]

也可以求解多变量的方程(组)：

  1. sage: x, y = var('x, y')
  2. sage: solve([x+y==6, x-y==4], x, y)
  3. [[x == 5, y == 1]]

下面是一个由Jason Grout提供的Sage求解非线性方程组的例子。我们先求方程组的符号解。

  1. sage: var('x y p q')
  2. (x, y, p, q)
  3. sage: eq1 = p+q==9
  4. sage: eq2 = q*y+p*x==-6
  5. sage: eq3 = q*y^2+p*x^2==24
  6. sage: solve([eq1,eq2,eq3,p==1],p,q,x,y)
  7. [[p == 1, q == 8, x == -4/3*sqrt(10) - 2/3, y == 1/6*sqrt(10) - 2/3],[p == 1, q == 8, x == 4/3*sqrt(10) - 2/3, y == -1/6*sqrt(10) - 2/3]]

要求解的近似值,可以这样：

  1. sage: solns = solve([eq1,eq2,eq3,p==1],p,q,x,y, solution_dict=True)
  2. sage:[[s[p].n(30), s[q].n(30), s[x].n(30), s[y].n(30)] for s in solns]
  3. [[1.0000000, 8.0000000, -4.8830369, -0.13962039],
  4.  [1.0000000, 8.0000000, 3.5497035, -1.1937129]]

(函数n输出数值的近似值,参数是以bit为单位的结果精度)

求方程的数值解

很多时候,solve找不到给定方程或方程组的精确解。这时候,你可以用find_root去找一个数值解。比如对于下面的方程,solve不会返回任何有用的信息：

  1. sage: theta = var('theta')
  2. sage: solve(cos(theta)==sin(theta), theta)
  3. [sin(theta) == cos(theta)]

但是我们可以用find_root在区间$0 < \phi < \pi/2$上寻找上述方程的解。

  1. sage: phi = var('phi')
  2. sage: find_root(cos(phi)==sin(phi),0,pi/2)
  3. 0.785398163397448...

微分、积分等

Sage知道如何求很多函数的微分和积分。比如求$\sin(u)$对$u$的微分这样做：

  1. sage: u = var('u')
  2. sage: diff(sin(u), u)
  3. cos(u)

计算$\sin(x^2)$的4阶微分：

  1. sage: diff(sin(x^2), x, 4)
  2. 16*x^4*sin(x^2) - 48*x^2*cos(x^2) - 12*sin(x^2)

分别计算$x^2+17y^2$对 x 和 y 的偏微分：

  1. sage: x, y = var('x,y')
  2. sage: f = x^2 + 17*y^2
  3. sage: f.diff(x)
  4. 2*x
  5. sage: f.diff(y)
  6. 34*y

再来看积分,定积分、不定积分都可以计算。来计算$\int x\sin(x^2)\, dx$和$\int_0^1 \frac{x}{x^2+1}\, dx$

  1. sage: integral(x*sin(x^2), x)
  2. -1/2*cos(x^2)
  3. sage: integral(x/(x^2+1), x, 0, 1)
  4. 1/2*log(2)

计算$\frac{1}{x^2-1}$的部分分式分解：

  1. sage: f = 1/((1+x)*(x-1))
  2. sage: f.partial_fraction(x)
  3. 1/2/(x - 1) - 1/2/(x + 1)

解微分方程组

可以用Sage求解常微分方程组。求解方程$x’+x-1=0$:

  1. sage: t = var('t')    # 定义变量t
  2. sage: x = function('x',t)   # 定义x是变量t的函数
  3. sage: DE = diff(x, t) + x - 1
  4. sage: desolve(DE,[x,t])
  5. (_C + e^t)*e^(-t)

在这里其实Sage调用了Maxima的接口,所以它的输出看上去与其他Sage的输出略有不同。这里,上述微分方程的通解是：$x(t) = e^{-t}(e^{t}+c)$.

你也可以计算Laplace变换。下面计算$t^2e^t -\sin(t)$的Laplace变换：

  1. sage: s = var("s")
  2. sage: t = var("t")
  3. sage: f = t^2*exp(t) - sin(t)
  4. sage: f.laplace(t,s)
  5. -1/(s^2 + 1) + 2/(s - 1)^3

这儿有一个更复杂的例子。两个弹簧连在左边的墙上,

  1. |------\/\/\/\/\---|mass1|----\/\/\/\/\/----|mass2|
  2.          spring1               spring2

物体偏离平衡态的位移可以描述为一个2阶常微分方程：

m_1 x_1’’ + (k_1+k_2) x_1 - k_2 x_2 = 0 m_2 x_2’’+ k_2 (x_2-x_1) = 0,

这里$m{i}$是物体 i 的质量,$x{i}$是物体 i 偏离平衡态的位移,$k_{i}$是弹簧 i 的弹性系数。

例： 在下面的条件下,使用Sage求解上面的问题$m{1}=2$,$m{2}=1$,$k{1}=4$,$k{2}=2$,$x{1}(0)=3$,$x{1}’(0)=0$,$x{2}(0)=3$,$x{2}’(0)=0$.

解：对第一个方程做Laplace变换(记$x=x{1}$,$y=x{2}$)：

  1. sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
  2. sage: lde1 = de1.laplace("t","s"); lde1
  3. 2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)

结果很难读,意思其实是：

-2x’(0) + 2s^2*X(s) - 2sx(0) - 2Y(s) + 6X(s) = 0

(这里对函数$x(t)$的Laplace变换记为$X(s)$)。对第二个方程做Laplace变换：

  1. sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
  2. sage: lde2 = de2.laplace("t","s"); lde2
  3. (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s

即

-Y’(0) + s^2Y(s) + 2Y(s) - 2X(s) - sy(0) = 0

代入$x(0)$,$x’(0)$,$y(0)$, 和$y’(0)$的初始条件,并求解求出的两个方程：

  1. sage: var('s X Y')
  2. (s, X, Y)
  3. sage: eqns =[(2*s^2+6)*X-2*Y == 6*s, -2*X +(s^2+2)*Y == 3*s]
  4. sage: solve(eqns, X,Y)
  5. [[X == 3*(s^3 + 3*s)/(s^4 + 5*s^2 + 4), 
  6.   Y == 3*(s^3 + 5*s)/(s^4 + 5*s^2 + 4)]]

现在做逆Laplace变换得到结果：

  1. sage: var('s t')
  2. (s, t)
  3. sage: inverse_laplace((3*s^3 + 9*s)/(s^4 + 5*s^2 + 4),s,t)
  4. cos(2*t) + 2*cos(t)
  5. sage: inverse_laplace((3*s^3 + 15*s)/(s^4 + 5*s^2 + 4),s,t)
  6. -cos(2*t) + 4*cos(t)

所以,原方程组的解是：

x_1(t) = \cos(2t) + 2\cos(t), \quad x_2(t) = 4\cos(t) - \cos(2t).

可以把结果画出来：

  1. sage: t = var('t')
  2. sage: P = parametric_plot((cos(2*t) + 2*cos(t), 4*cos(t) - cos(2*t) ),\
  3. ...   (t, 0, 2*pi), rgbcolor=hue(0.9))
  4. sage: show(P)

每一个分支都可以画出来：

  1. sage: t = var('t')
  2. sage: p1 = plot(cos(2*t) + 2*cos(t), (t,0, 2*pi), rgbcolor=hue(0.3))
  3. sage: p2 = plot(4*cos(t) - cos(2*t), (t,0, 2*pi), rgbcolor=hue(0.6))
  4. sage: show(p1 + p2)

更多关于做图的内容,参见绘图.

可以阅读[NagleEtAl2004][]的第5.5节来了解微分方程。

解微分方程组的Euler方法

下面的例子中,我们展示求解1阶,2阶常微分方程组的Euler方法。我们先来回顾一下1阶方程的基本知识。给定如下形式的初值问题：

y’=f(x,y),y(a)=c

我们要找方程在$x=b$处的近似解,且$b>a$.

根据微分的定义

y’(x) \approx \frac{y(x+h)-y(x)}{h},

这里$h>0$是给定的,且较小的量。与微分方程一起得到$f(x,y(x))\approx\frac{y(x+h)-y(x)}{h}$. 现在求$y(x+h)$:

y(x+h) \approx y(x) + h*f(x,y(x)).

如果将$h f(x,y(x))$称为”校正项”(没有更好的名字), 称$y(x)$为$y$的旧值,$y(x+h)$为$y$的新值, 那么该近似公式可以改写为：

y{new} \approx y{old} + h*f(x,y_{old}).

如果将由a到b的区间n等分,则$h=\frac{b-a}{n}$,我们可以用一个表记录该方法得到的信息。

$x$ $y$ $hf(x,y)$
$a$ $c$ $hf(a,c)$
$a+h$ $c+hf(a,c)$ …
$a+2h$ …
…
$b=a+nh$ ??? …

我们的目标是把表中的空格都填上,每次一行,直到到达 ??? 这一项, 也就是Euler方法求得的$y(b)$的近似值。

类似的方法可以来求解常微分方程组。

例：用4步Euler方法求$z(t)$在$t=1$处的近似值, 这里$z’’+tz’+z=0$,$z(0)=1$,$z’(0)=0$.

我们必须将2阶常微分方程化为两个1阶微分方程(令$x=z$,$y=z’$)并再应用Euler方法：

  1. sage: t,x,y = PolynomialRing(RealField(10),3,"txy").gens()
  2. sage: f = y; g = -x - y * t
  3. sage: eulers_method_2x2(f,g, 0, 1, 0, 1/4, 1)
  4.       t                x            h*f(t,x,y)                y       h*g(t,x,y)
  5.       0                1                  0.00                0           -0.25
  6.     1/4              1.0                -0.062            -0.25           -0.23
  7.     1/2             0.94                 -0.12            -0.48           -0.17
  8.     3/4             0.82                 -0.16            -0.66          -0.081
  9.       1             0.65                 -0.18            -0.74           0.022

即,$z(1)\approx 0.65$.

我们可以把点$(x,y)$画出来,得到曲线的近似图像。 函数eulers_method_2x2_plot可以做到这一点。 为了应用该函数,要先定义函数$f$和$g$来接受含三个坐标的参数： ($t, x, y$).

  1. sage: f = lambda z: z[2]       # f(t,x,y) = y
  2. sage: g = lambda z: -sin(z[1])  # g(t,x,y) = -sin(x)
  3. sage: P = eulers_method_2x2_plot(f,g, 0.0, 0.75, 0.0, 0.1, 1.0)

这里,P保存了两个图像,P[0]是 $x$ 关于 $t$ 的图像,P[1]是 $y$关于 $t$ 的图像。我们把它们都画出来：

  1. sage: show(P[0]+ P[1])

(更多关于做图的内容,参见绘图.)

特殊函数

一些正交多项式和特殊函数是使用PARIGAP 和MaximaMax 实现的。在Sage参考手册的相关章节(“正交多项式”和”特殊函数”)中有详细信息。

  1. sage: x = polygen(QQ, 'x')
  2. sage: chebyshev_U(2,x)
  3. 4*x^2 - 1
  4. sage: bessel_I(1,1,"pari",250)
  5. 0.56515910399248502720769602760986330732889962162109200948029448947925564096
  6. sage: bessel_I(1,1).n()
  7. 0.565159103992485
  8. sage: bessel_I(2,1.1).n
  9. 0.167089499251049

这里Sage直接求得数值解,如果想求符号解,请像下面这样直接使用Maxima接口：

  1. sage: maxima.eval("f:bessel_y(v, w)")
  2. 'bessel_y(v,w)'
  3. sage: maxima.eval("diff(f,w)")
  4. '(bessel_y(v-1,w)-bessel_y(v+1,w))/2'

向量

参见Vector Calculus Tutorial

[NagleEtAl2004]:Nagle, Saff, and Snider. Fundamentals of Differential Equations. 6th edition, Addison-Wesley,2004.