基础代数与微积分

SAGE可以执行与基本代数和微积分相关的各种计算:例如,寻找方程的解、微分、积分和拉普拉斯变换。请参阅 Sage Constructions 有关更多示例的文档。

在所有这些示例中,重要的是要注意,函数中的变量定义为 var(...) 。举个例子:

sage: u = var('u')
sage: diff(sin(u), u)
cos(u)

如果你得到了一个 NameError ,检查您是否拼错了什么,或者忘记用定义变量 var(...)

解方程组

精确求解方程

这个 solve 函数求解方程。要使用它,首先指定一些变量;然后将参数 solve 是一个方程(或一组方程),以及要求解的变量:

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

您可以用其他变量来求解一个变量的方程式:

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

您还可以求解以下几个变量:

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

Jason Grout提供了使用Sage求解非线性方程组的以下示例:首先,我们对该系统进行符号求解:

sage: var('x y p q')
(x, y, p, q)
sage: eq1 = p+q==9
sage: eq2 = q*y+p*x==-6
sage: eq3 = q*y^2+p*x^2==24
sage: solve([eq1,eq2,eq3,p==1],p,q,x,y)
[[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]]

对于解的数值近似,您可以改用:

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

(该功能 n 打印数值近似值,参数是精度位数。)

用数值方法求解方程

很多时候, solve 将无法找到指定的一个或多个方程的精确解。当失败时,您可以使用 find_root 找到一个数值解。例如,Solve不会为以下方程式返回任何有趣的内容:

sage: theta = var('theta')
sage: solve(cos(theta)==sin(theta), theta)
[sin(theta) == cos(theta)]

另一方面,我们可以使用 find_root 要在这个范围内找到上述方程的解 \(0 < \phi < \pi/2\) **

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

差异化、一体化等

Sage知道如何区分和集成许多功能。例如,为了区分 \(\sin(u)\) 关于… \(u\) ,请执行以下操作:

sage: u = var('u')
sage: diff(sin(u), u)
cos(u)

计算四阶导数 \(\sin(x^2)\)

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

计算的偏导数 \(x^2+17y^2\) 关于… xy ,分别为:

sage: x, y = var('x,y')
sage: f = x^2 + 17*y^2
sage: f.diff(x)
2*x
sage: f.diff(y)
34*y

We move on to integrals, both indefinite and definite. To compute \(\int x\sin(x^2)\, dx\) and \(\int_0^1 \frac{x}{x^2+1}\, dx\)

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

计算部分分数分解的步骤 \(\frac{1}{x^2-1}\)

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

解微分方程式

您可以使用Sage研究常微分方程式。要解这个方程 \(x'+x-1=0\)

sage: t = var('t')    # define a variable t
sage: x = function('x')(t)   # define x to be a function of that variable
sage: DE = diff(x, t) + x - 1
sage: desolve(DE, [x,t])
(_C + e^t)*e^(-t)

这使用Sage到Maxima的接口 [Max], 因此,它的输出可能与其他Sage输出略有不同。在这种情况下,这意味着微分方程式的一般解是 \(x(t) = e^{-t}(e^{t}+c)\)

你也可以计算拉普拉斯变换;拉普拉斯变换 \(t^2e^t -\sin(t)\) 按如下方式计算:

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

这里有一个更复杂的例子。连接在左侧墙上的耦合弹簧偏离平衡的位移

|------\/\/\/\/\---|mass1|----\/\/\/\/\/----|mass2|
         spring1               spring2

由二阶微分方程组模拟

\[ \begin{align}\begin{aligned}M_1 x_1‘’+(k_1+k_2)x_1-k_2 x_2=0\\M_2x_2‘’+k_2(x_2-x_1)=0,\end{aligned}\end{align} \]

哪里 \(m_{i}\) 是物体的质量 i\(x_{i}\) 是从质量平衡的位移 i ,以及 \(k_{i}\) 春天是春天的常量吗? i

Example: 用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\)

解:取第一个方程的拉普拉斯变换(带符号 \(x=x_{1}\)\(y=x_{2}\) ):

sage: t,s = SR.var('t,s')
sage: x = function('x')
sage: y = function('y')
sage: f = 2*x(t).diff(t,2) + 6*x(t) - 2*y(t)
sage: f.laplace(t,s)
2*s^2*laplace(x(t), t, s) - 2*s*x(0) + 6*laplace(x(t), t, s) - 2*laplace(y(t), t, s) - 2*D[0](x)(0)

这很难读懂,但上面写着

\[-2x‘(0)+2s^2\cot X(S)-2sx(0)-2Y(S)+6X(S)=0\]

(其中,小写字母的拉普拉斯变换如下 \(x(t)\) 是大写字母函数 \(X(s)\) )。以第二个方程的拉普拉斯变换为例:

sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
sage: lde2 = de2.laplace("t","s"); lde2.sage()
s^2*laplace(y(t), t, s) - s*y(0) - 2*laplace(x(t), t, s) + 2*laplace(y(t), t, s) - D[0](y)(0)

这上面写着

\[-Y‘(0)+S^2Y(S)+2Y(S)-2X(S)-SY(0)=0。\]

插入以下初始条件: \(x(0)\)\(x'(0)\)\(y(0)\) ,以及 \(y'(0)\) ,并求解得到的两个方程:

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

现在用拉普拉斯逆变换来得到答案:

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

因此,解决方案是

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

可以使用以下命令以参数方式绘制该图

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

可以使用以下工具绘制各个组件

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

有关打印的详细信息,请参见 标绘 。请参阅第5.5节 [NagleEtAl2004] 以获取有关微分方程的更多信息。

微分方程组的欧拉法

在下一个例子中,我们将说明一阶和二阶常微分方程组的欧拉方法。我们首先回顾一下一阶方程的基本思想。给定一个形式的初值问题

\[Y‘=f(x,y),\quad y(A)=c,\]

我们希望在以下位置找到解的近似值 \(x=b\) 使用 \(b>a\)

从衍生品的定义中回忆起

\[Y‘(X)\约x\frc{y(x+h)-y(X)}{h},\]

哪里 \(h>0\) 是给与的,而且很小。这一点和DE一起提供了 \(f(x,y(x))\approx \frac{y(x+h)-y(x)}{h}\) 。现在解出 \(y(x+h)\)

\[Y(x+h)\近似y(X)+h\cot f(x,y(X))。\]

如果我们打电话给 \(h \cdot f(x,y(x))\) “修正项”(因为没有更好的东西),调用 \(y(x)\) “旧价值” y “,并呼叫 \(y(x+h)\) 的“新价值” y 则该近似值可以重新表示为

\[Y_{new}\近似y_{old}+h\cot f(x,y_{old})。\]

如果我们将间隔从 ab vt.进入,进入 n 步骤,这样就可以 \(h=\frac{b-a}{n}\) ,然后我们可以将此方法的信息记录在表中。

\(x\)

\(y\)

\(h\cdot f(x,y)\)

\(a\)

\(c\)

\(h\cdot f(a,c)\)

\(a+h\)

\(c+h\cdot f(a,c)\)

..。

\(a+2h\)

..。

..。

\(b=a+nh\)

???

..。

我们的目标是填充表中的所有空白处,一次一行,直到我们达到?Entry,这是Euler方法对 \(y(b)\)

对于ODE系统的想法也是类似的。

Example: 数值近似值 \(z(t)\) 在… \(t=1\) 使用欧拉法的4个步骤,其中 \(z''+tz'+z=0\)\(z(0)=1\)\(z'(0)=0\)

我们必须将二阶常数降为两个一阶DES的系统(使用 \(x=z\)\(y=z'\) ),并应用欧拉方法:

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

因此, \(z(1)\approx 0.65\)

我们也可以画出这些点 \((x,y)\) 以获得曲线的大致图像。该功能 eulers_method_2x2_plot 为了使用它,我们需要定义函数 fg 它接受一个带有三个坐标的参数: (txy )。

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

在这点上, P 正在存储两个绘图: P[0] ,故事的情节 x v.v. t ,以及 P[1] ,故事的情节 y v.v. t 。我们可以将这两种情况绘制如下:

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

(有关绘制的更多信息,请参见 标绘 。)

特殊功能

用这两个参数实现了几个正交多项式和特殊函数 [GAP] 和Maxima [Max]. 在Sage参考手册的相应章节(分别为“正交多项式”和“特殊函数”)中记录了这些内容。

sage: x = polygen(QQ, 'x')
sage: chebyshev_U(2,x)
4*x^2 - 1
sage: bessel_I(1,1).n(250)
0.56515910399248502720769602760986330732889962162109200948029448947925564096
sage: bessel_I(1,1).n()
0.565159103992485
sage: bessel_I(2,1.1).n()
0.167089499251049

在这一点上,Sage仅包装了这些函数用于数值用途。对于符号用途,请直接使用Maxima接口,如下例所示:

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

向量微积分

请参阅 Vector Calculus Tutorial