基础代数与微积分¶
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]]
以下使用Sage求解非线性方程组的示例由Jason screet提供:首先,我们象征性地求解该系统:
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\) 关于 x 和 y ,分别为:
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
我们继续学习积分,包括不定积分和定积分。计算 \(\int x\sin(x^2)\, dx\) 和 \(\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
是由二阶微分方程组建立的
在哪里? \(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\) .
解:取第一个方程的拉普拉斯变换(带符号 \(x=x_{{1}}\) , \(y=x_{{2}}\) ):
sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
sage: lde1 = de1.laplace("t","s"); lde1
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)
这很难读懂,但上面写着
(其中小写函数的拉普拉斯变换 \(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
(-%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
上面写着
插入初始条件 \(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)
因此,解决办法是
这可以用参数化绘图
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] 关于微分方程的更多信息。
微分方程组的欧拉方法¶
在下一个例子中,我们将说明一阶和二阶常微分方程的欧拉方法。我们首先回顾一下一阶方程的基本思想。给出了一个形式的初值问题
我们想在 \(x=b\) 具有 \(b>a\) .
回想一下导数的定义
在哪里? \(h>0\) 是给予和小。这个和DE一起给 \(f(x,y(x))\approx \frac{{y(x+h)-y(x)}}{{h}}\) . 现在解决 \(y(x+h)\) :
如果我们打电话 \(h \cdot f(x,y(x))\) “修正条款”(因为没有更好的东西),调用 \(y(x)\) “旧价值” y “,然后打电话 \(y(x+h)\) “新价值 y ,则该近似值可重新表示为
如果我们从 a 到 b 进入之内 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\) |
???? |
… |
我们的目标是填写表格中的所有空格,一次一行,直到我们到达???入口,这是Euler方法的近似值 \(y(b)\) .
《诗经》的理念与此相似。
例子: 数值近似 \(z(t)\) 在 \(t=1\) 使用欧拉方法的4个步骤,其中 \(z''+tz'+z=0\) , \(z(0)=1\) , \(z'(0)=0\) .
我们必须将二阶常微分方程降为两个一阶微分方程组(使用 \(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
将这样做;为了使用它,我们需要定义函数 f 和 g 一个参数有三个坐标: (t , x , y )
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 VS t 和 P[1]
,情节 y VS t . 我们可以把这两个图都画出来:
sage: show(P[0] + P[1])
(有关打印的详细信息,请参见 作图 )
特殊功能¶
利用这两种方法实现了几个正交多项式和特殊函数 [GAP] 还有马克西玛 [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'