|
| 1 | +******** |
| 2 | +微积分 |
| 3 | +******** |
| 4 | + |
| 5 | +这里有一些使用 Sage 进行微积分符号计算的例子。 |
| 6 | + |
| 7 | + |
| 8 | +.. index:: |
| 9 | + pair: calculus; differentiation |
| 10 | + |
| 11 | +微分 |
| 12 | +=============== |
| 13 | + |
| 14 | +微分: |
| 15 | + |
| 16 | +:: |
| 17 | + |
| 18 | + sage: var('x k w') |
| 19 | + (x, k, w) |
| 20 | + sage: f = x^3 * e^(k*x) * sin(w*x); f |
| 21 | + x^3*e^(k*x)*sin(w*x) |
| 22 | + sage: f.diff(x) |
| 23 | + w*x^3*cos(w*x)*e^(k*x) + k*x^3*e^(k*x)*sin(w*x) + 3*x^2*e^(k*x)*sin(w*x) |
| 24 | + sage: latex(f.diff(x)) |
| 25 | + w x^{3} \cos\left(w x\right) e^{\left(k x\right)} + k x^{3} e^{\left(k x\right)} \sin\left(w x\right) + 3 \, x^{2} e^{\left(k x\right)} \sin\left(w x\right) |
| 26 | + |
| 27 | +如果你输入 ``view(f.diff(x))``,将会打开一个新的窗口显示编译后的输出。在 notebook 中,你可以在单元格中输入 |
| 28 | + |
| 29 | +.. CODE-BLOCK:: ipython |
| 30 | +
|
| 31 | + var('x k w') |
| 32 | + f = x^3 * e^(k*x) * sin(w*x) |
| 33 | + show(f) |
| 34 | + show(f.diff(x)) |
| 35 | +
|
| 36 | +然后按 ``shift-enter`` 来获得类似的结果。 |
| 37 | + |
| 38 | +你还可以在 notebook 单元格中,使用以下命令进行微分和积分 |
| 39 | + |
| 40 | +.. CODE-BLOCK:: ipython |
| 41 | +
|
| 42 | + R = PolynomialRing(QQ,"x") |
| 43 | + x = R.gen() |
| 44 | + p = x^2 + 1 |
| 45 | + show(p.derivative()) |
| 46 | + show(p.integral()) |
| 47 | +
|
| 48 | +或者在命令行中,使用以下命令进行微分和积分 |
| 49 | + |
| 50 | +:: |
| 51 | + |
| 52 | + sage: R = PolynomialRing(QQ,"x") |
| 53 | + sage: x = R.gen() |
| 54 | + sage: p = x^2 + 1 |
| 55 | + sage: p.derivative() |
| 56 | + 2*x |
| 57 | + sage: p.integral() |
| 58 | + 1/3*x^3 + x |
| 59 | + |
| 60 | +此时你也可以输入 ``view(p.derivative())`` 或 ``view(p.integral())`` 以打开一个新窗口, |
| 61 | +并由 LaTeX 对输出进行排版。 |
| 62 | + |
| 63 | +.. index:: |
| 64 | + pair: calculus; critical points |
| 65 | + |
| 66 | +临界点 |
| 67 | +--------------- |
| 68 | + |
| 69 | +你可以找到分段定义函数的临界点: |
| 70 | + |
| 71 | +:: |
| 72 | + |
| 73 | + sage: x = PolynomialRing(RationalField(), 'x').gen() |
| 74 | + sage: f1 = x^0 |
| 75 | + sage: f2 = 1-x |
| 76 | + sage: f3 = 2*x |
| 77 | + sage: f4 = 10*x-x^2 |
| 78 | + sage: f = piecewise([((0,1),f1), ((1,2),f2), ((2,3),f3), ((3,10),f4)]) |
| 79 | + sage: f.critical_points() |
| 80 | + [5.0] |
| 81 | + |
| 82 | +.. index:: Taylor series, power series |
| 83 | + |
| 84 | +幂级数 |
| 85 | +------------ |
| 86 | + |
| 87 | +Sage 提供了几种构建和处理幂级数的方法。 |
| 88 | + |
| 89 | +为了从函数表达式中获取泰勒级数,请在表达式上使用 ``.taylor()`` 方法:: |
| 90 | + |
| 91 | + sage: var('f0 k x') |
| 92 | + (f0, k, x) |
| 93 | + sage: g = f0/sinh(k*x)^4 |
| 94 | + sage: g.taylor(x, 0, 3) |
| 95 | + -62/945*f0*k^2*x^2 + 11/45*f0 - 2/3*f0/(k^2*x^2) + f0/(k^4*x^4) |
| 96 | + |
| 97 | +可以通过 ``.series()`` 方法获得函数的形式幂级数展开:: |
| 98 | + |
| 99 | + sage: (1/(2-cos(x))).series(x,7) |
| 100 | + 1 + (-1/2)*x^2 + 7/24*x^4 + (-121/720)*x^6 + Order(x^7) |
| 101 | + |
| 102 | +然而,目前对此类序列的某些操作很难执行。有两种替代方案:要么使用 Sage 的 Maxima 子系统以获取完整符号功能:: |
| 103 | + |
| 104 | + sage: f = log(sin(x)/x) |
| 105 | + sage: f.taylor(x, 0, 10) |
| 106 | + -1/467775*x^10 - 1/37800*x^8 - 1/2835*x^6 - 1/180*x^4 - 1/6*x^2 |
| 107 | + sage: maxima(f).powerseries(x,0)._sage_() |
| 108 | + sum(2^(2*i... - 1)*(-1)^i...*x^(2*i...)*bern(2*i...)/(i...*factorial(2*i...)), i..., 1, +Infinity) |
| 109 | + |
| 110 | +要么你可以使用形式幂级数环进行快速计算。这些环缺乏符号函数:: |
| 111 | + |
| 112 | + sage: R.<w> = QQ[[]] |
| 113 | + sage: ps = w + 17/2*w^2 + 15/4*w^4 + O(w^6); ps |
| 114 | + w + 17/2*w^2 + 15/4*w^4 + O(w^6) |
| 115 | + sage: ps.exp() |
| 116 | + 1 + w + 9*w^2 + 26/3*w^3 + 265/6*w^4 + 413/10*w^5 + O(w^6) |
| 117 | + sage: (1+ps).log() |
| 118 | + w + 8*w^2 - 49/6*w^3 - 193/8*w^4 + 301/5*w^5 + O(w^6) |
| 119 | + sage: (ps^1000).coefficients() |
| 120 | + [1, 8500, 36088875, 102047312625, 1729600092867375/8] |
| 121 | + |
| 122 | +.. index:: |
| 123 | + pair: calculus; integration |
| 124 | + |
| 125 | +积分 |
| 126 | +=========== |
| 127 | + |
| 128 | +下面的 :ref:`section-riemannsums` 讨论了数值积分。 |
| 129 | + |
| 130 | +Sage 可以自行对一些简单函数进行积分: |
| 131 | + |
| 132 | +:: |
| 133 | + |
| 134 | + sage: f = x^3 |
| 135 | + sage: f.integral(x) |
| 136 | + 1/4*x^4 |
| 137 | + sage: integral(x^3,x) |
| 138 | + 1/4*x^4 |
| 139 | + sage: f = x*sin(x^2) |
| 140 | + sage: integral(f,x) |
| 141 | + -1/2*cos(x^2) |
| 142 | + |
| 143 | +Sage 还可以计算涉及极限的符号定积分。 |
| 144 | + |
| 145 | +:: |
| 146 | + |
| 147 | + sage: var('x, k, w') |
| 148 | + (x, k, w) |
| 149 | + sage: f = x^3 * e^(k*x) * sin(w*x) |
| 150 | + sage: f.integrate(x) |
| 151 | + ((24*k^3*w - 24*k*w^3 - (k^6*w + 3*k^4*w^3 + 3*k^2*w^5 + w^7)*x^3 + 6*(k^5*w + 2*k^3*w^3 + k*w^5)*x^2 - 6*(3*k^4*w + 2*k^2*w^3 - w^5)*x)*cos(w*x)*e^(k*x) - (6*k^4 - 36*k^2*w^2 + 6*w^4 - (k^7 + 3*k^5*w^2 + 3*k^3*w^4 + k*w^6)*x^3 + 3*(k^6 + k^4*w^2 - k^2*w^4 - w^6)*x^2 - 6*(k^5 - 2*k^3*w^2 - 3*k*w^4)*x)*e^(k*x)*sin(w*x))/(k^8 + 4*k^6*w^2 + 6*k^4*w^4 + 4*k^2*w^6 + w^8) |
| 152 | + sage: integrate(1/x^2, x, 1, infinity) |
| 153 | + 1 |
| 154 | + |
| 155 | + |
| 156 | +.. index: convolution |
| 157 | +
|
| 158 | +卷积 |
| 159 | +----------- |
| 160 | + |
| 161 | +你可以计算任意分段函数与另一个函数的卷积(在定义域之外,它们被假定为零)。 |
| 162 | +以下是 :math:`f`, :math:`f*f` 和 :math:`f*f*f` 的定义, |
| 163 | +其中 :math:`f(x)=1`, :math:`0<x<1`: |
| 164 | + |
| 165 | +:: |
| 166 | + |
| 167 | + sage: x = PolynomialRing(QQ, 'x').gen() |
| 168 | + sage: f = piecewise([((0,1),1*x^0)]) |
| 169 | + sage: g = f.convolution(f) |
| 170 | + sage: h = f.convolution(g) |
| 171 | + sage: set_verbose(-1) |
| 172 | + sage: P = f.plot(); Q = g.plot(rgbcolor=(1,1,0)); R = h.plot(rgbcolor=(0,1,1)) |
| 173 | + |
| 174 | +要查看此内容,请输入 ``show(P+Q+R)``。 |
| 175 | + |
| 176 | + |
| 177 | +.. _section-riemannsums: |
| 178 | + |
| 179 | +黎曼和与梯形法积分 |
| 180 | +---------------------------------------- |
| 181 | + |
| 182 | +关于 :math:`\int_a^bf(x)\, dx` 的数值近似, |
| 183 | +其中 :math:`f` 是分段函数,可以: |
| 184 | + |
| 185 | + |
| 186 | +- 计算(用于绘图目的)根据梯形法则定义的分段线性函数,基于将其分割为 :math:`N` 个子区间进行数值积分; |
| 187 | + |
| 188 | +- 梯形法则给出的近似值; |
| 189 | + |
| 190 | +- 计算(用于绘图目的)根据黎曼和(左端点、右端点或中点)定义的分段常数函数, |
| 191 | + 基于将其分割为 :math:`N` 个子区间进行数值积分; |
| 192 | + |
| 193 | +- 黎曼和近似值给出的近似值。 |
| 194 | + |
| 195 | + |
| 196 | +:: |
| 197 | + |
| 198 | + sage: f1(x) = x^2 |
| 199 | + sage: f2(x) = 5-x^2 |
| 200 | + sage: f = piecewise([[[0,1], f1], [RealSet.open_closed(1,2), f2]]) |
| 201 | + sage: t = f.trapezoid(2); t |
| 202 | + piecewise(x|-->1/2*x on (0, 1/2), x|-->3/2*x - 1/2 on (1/2, 1), x|-->7/2*x - 5/2 on (1, 3/2), x|-->-7/2*x + 8 on (3/2, 2); x) |
| 203 | + sage: t.integral() |
| 204 | + piecewise(x|-->1/4*x^2 on (0, 1/2), x|-->3/4*x^2 - 1/2*x + 1/8 on (1/2, 1), x|-->7/4*x^2 - 5/2*x + 9/8 on (1, 3/2), x|-->-7/4*x^2 + 8*x - 27/4 on (3/2, 2); x) |
| 205 | + sage: t.integral(definite=True) |
| 206 | + 9/4 |
| 207 | + |
| 208 | +.. index: Laplace transform |
| 209 | +
|
| 210 | +拉普拉斯变换 |
| 211 | +------------------ |
| 212 | + |
| 213 | +如果你有一个分段定义的多项式函数,那么有一个“原生”命令用于计算拉普拉斯变换。 |
| 214 | +这将调用 Maxima,但值得注意的是,Maxima 无法(使用最后几个示例中的直接接口)处理这种类型的计算。 |
| 215 | + |
| 216 | +:: |
| 217 | + |
| 218 | + sage: var('x s') |
| 219 | + (x, s) |
| 220 | + sage: f1(x) = 1 |
| 221 | + sage: f2(x) = 1-x |
| 222 | + sage: f = piecewise([((0,1),f1), ((1,2),f2)]) |
| 223 | + sage: f.laplace(x, s) |
| 224 | + -e^(-s)/s + (s + 1)*e^(-2*s)/s^2 + 1/s - e^(-s)/s^2 |
| 225 | + |
| 226 | +对于其他“合理”的函数,可以使用 Maxima 接口计算拉普拉斯变换: |
| 227 | + |
| 228 | +:: |
| 229 | + |
| 230 | + sage: var('k, s, t') |
| 231 | + (k, s, t) |
| 232 | + sage: f = 1/exp(k*t) |
| 233 | + sage: f.laplace(t,s) |
| 234 | + 1/(k + s) |
| 235 | + |
| 236 | +上面是计算拉普拉斯变换的一种方法 |
| 237 | + |
| 238 | +:: |
| 239 | + |
| 240 | + sage: var('s, t') |
| 241 | + (s, t) |
| 242 | + sage: f = t^5*exp(t)*sin(t) |
| 243 | + sage: L = laplace(f, t, s); L |
| 244 | + 3840*(s - 1)^5/(s^2 - 2*s + 2)^6 - 3840*(s - 1)^3/(s^2 - 2*s + 2)^5 + |
| 245 | + 720*(s - 1)/(s^2 - 2*s + 2)^4 |
| 246 | + |
| 247 | +上面是另一种方法。 |
| 248 | + |
| 249 | +.. index: |
| 250 | + pair: differential equations; solve |
| 251 | +
|
| 252 | +常微分方程 |
| 253 | +=============================== |
| 254 | + |
| 255 | +使用 Sage 接口与 Maxima 可以符号化地求解常微分方程。参见 |
| 256 | + |
| 257 | +:: |
| 258 | + |
| 259 | + sage:desolvers? |
| 260 | + |
| 261 | +获取可用命令。 |
| 262 | +可以使用 Sage 接口与 Octave(一个实验性包)或 GSL(Gnu 科学库)中的例程来数值求解常微分方程。 |
| 263 | + |
| 264 | +例如,通过 Sage 的 Maxima 接口符号化地求解常微分方程(请勿输入 ``....:``): |
| 265 | + |
| 266 | +:: |
| 267 | + |
| 268 | + sage: y=function('y')(x); desolve(diff(y,x,2) + 3*x == y, dvar = y, ics = [1,1,1]) |
| 269 | + 3*x - 2*e^(x - 1) |
| 270 | + sage: desolve(diff(y,x,2) + 3*x == y, dvar = y) |
| 271 | + _K2*e^(-x) + _K1*e^x + 3*x |
| 272 | + sage: desolve(diff(y,x) + 3*x == y, dvar = y) |
| 273 | + (3*(x + 1)*e^(-x) + _C)*e^x |
| 274 | + sage: desolve(diff(y,x) + 3*x == y, dvar = y, ics = [1,1]).expand() |
| 275 | + 3*x - 5*e^(x - 1) + 3 |
| 276 | + |
| 277 | + sage: f=function('f')(x); desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f, ics = [0,1,2]) |
| 278 | + x*e^x + e^x |
| 279 | + |
| 280 | + sage: desolve_laplace(diff(f,x,2) == 2*diff(f,x)-f, dvar = f) |
| 281 | + -x*e^x*f(0) + x*e^x*D[0](f)(0) + e^x*f(0) |
| 282 | + |
| 283 | +.. index: |
| 284 | + pair: differential equations; plot |
| 285 | +
|
| 286 | +如果你已经安装了 ``Octave`` 和 ``gnuplot``, |
| 287 | + |
| 288 | +:: |
| 289 | + |
| 290 | + sage: octave.de_system_plot(['x+y','x-y'], [1,-1], [0,2]) # optional - octave |
| 291 | + |
| 292 | +将在同一个图中绘制常微分方程组的两个图像 :math:`(t,x(t)), (t,y(t))` (:math:`t`-轴为横轴) |
| 293 | + |
| 294 | +.. math:: |
| 295 | + x' = x+y, x(0) = 1; y' = x-y, y(0) = -1, |
| 296 | +
|
| 297 | +对于 :math:`0 <= t <= 2`。使用 ``desolve_system_rk4`` 也可以获得相同的结果:: |
| 298 | + |
| 299 | + sage: x, y, t = var('x y t') |
| 300 | + sage: P=desolve_system_rk4([x+y, x-y], [x,y], ics=[0,1,-1], ivar=t, end_points=2) |
| 301 | + sage: p1 = list_plot([[i,j] for i,j,k in P], plotjoined=True) |
| 302 | + sage: p2 = list_plot([[i,k] for i,j,k in P], plotjoined=True, color='red') |
| 303 | + sage: p1+p2 |
| 304 | + Graphics object consisting of 2 graphics primitives |
| 305 | + |
| 306 | +该方程组也可以通过使用命令 ``desolve_system`` 来求解。 |
| 307 | + |
| 308 | +.. skip |
| 309 | +
|
| 310 | +:: |
| 311 | + |
| 312 | + sage: t=var('t'); x=function('x',t); y=function('y',t) |
| 313 | + sage: des = [diff(x,t) == x+y, diff(y,t) == x-y] |
| 314 | + sage: desolve_system(des, [x,y], ics = [0, 1, -1]) |
| 315 | + [x(t) == cosh(sqrt(2)*t), y(t) == sqrt(2)*sinh(sqrt(2)*t) - cosh(sqrt(2)*t)] |
| 316 | + |
| 317 | +此命令的输出 *不* 是一对函数。 |
| 318 | + |
| 319 | +最后,可以使用幂级数求解线性微分方程: |
| 320 | + |
| 321 | +:: |
| 322 | + |
| 323 | + sage: R.<t> = PowerSeriesRing(QQ, default_prec=10) |
| 324 | + sage: a = 2 - 3*t + 4*t^2 + O(t^10) |
| 325 | + sage: b = 3 - 4*t^2 + O(t^7) |
| 326 | + sage: f = a.solve_linear_de(prec=5, b=b, f0=3/5) |
| 327 | + sage: f |
| 328 | + 3/5 + 21/5*t + 33/10*t^2 - 38/15*t^3 + 11/24*t^4 + O(t^5) |
| 329 | + sage: f.derivative() - a*f - b |
| 330 | + O(t^4) |
| 331 | + |
| 332 | +周期函数的傅里叶级数 |
| 333 | +==================================== |
| 334 | + |
| 335 | +设 :math:`f` 是一个周期为 `2L` 的实周期函数。 |
| 336 | +`f` 的傅里叶级数是 |
| 337 | + |
| 338 | +.. MATH:: |
| 339 | +
|
| 340 | + S(x) = \frac{a_0}{2} + \sum_{n=1}^\infty \left[a_n\cos\left(\frac{n\pi x}{L}\right) + |
| 341 | + b_n\sin\left(\frac{n\pi x}{L}\right)\right] |
| 342 | +
|
| 343 | +其中 |
| 344 | + |
| 345 | +.. MATH:: |
| 346 | +
|
| 347 | + a_n = \frac{1}{L}\int_{-L}^L |
| 348 | + f(x)\cos\left(\frac{n\pi x}{L}\right) dx, |
| 349 | +
|
| 350 | +并且 |
| 351 | + |
| 352 | +.. MATH:: |
| 353 | +
|
| 354 | + b_n = \frac{1}{L}\int_{-L}^L |
| 355 | + f(x)\sin\left(\frac{n\pi x}{L}\right) dx, |
| 356 | +
|
| 357 | +傅里叶系数 `a_n` 和 `b_n` 是通过声明 `f` 在一个周期内分段定义的函数并调用方法 |
| 358 | +``fourier_series_cosine_coefficient`` 和 ``fourier_series_sine_coefficient`` 来计算的, |
| 359 | +而部分和是通过 ``fourier_series_partial_sum`` 获得的:: |
| 360 | + |
| 361 | + sage: f = piecewise([((0,pi/2), -1), ((pi/2,pi), 2)]) |
| 362 | + sage: f.fourier_series_cosine_coefficient(0) |
| 363 | + 1 |
| 364 | + sage: f.fourier_series_sine_coefficient(5) |
| 365 | + -6/5/pi |
| 366 | + sage: s5 = f.fourier_series_partial_sum(5); s5 |
| 367 | + -6/5*sin(10*x)/pi - 2*sin(6*x)/pi - 6*sin(2*x)/pi + 1/2 |
| 368 | + sage: plot(f, (0,pi)) + plot(s5, (x,0,pi), color='red') |
| 369 | + Graphics object consisting of 2 graphics primitives |
| 370 | + |
| 371 | +.. PLOT:: |
| 372 | + |
| 373 | + f = piecewise([((0,pi/2), -1), ((pi/2,pi), 2)]) |
| 374 | + s5 = f.fourier_series_partial_sum(5) |
| 375 | + g = plot(f, (0,pi)) + plot(s5, (x,0,pi), color='red') |
| 376 | + sphinx_plot(g) |
0 commit comments