# -*- coding:utf-8 -*-
import visual as vs
import ode

vs.scene.center = (0, 2, 0)

"""  ODE環境の設定  """
# ODEワールドの作成 + ワールドの重力設定
world = ode.World()
world.setGravity((0, -9.81, 0))

# 剛体の作成 + 質量特性の設定
Ball_Body = ode.Body(world) 
M = ode.Mass()              
M.setSphere(2500.0, 0.5)  
Ball_Body.setMass(M)        
Ball_Body.setPosition((0, 4, 0))

""" VPython環境の設定 """
Floor = vs.box(length=4, height=0.2, width=4)
Floor.pos = vs.vector(0, 0, 0)

Ball_VS = vs.sphere(pos=Ball_Body.getPosition(),
                       radius=0.5,
                       color=vs.color.red)

if __name__=="__main__":
    dt = 0.05
    total_time = 0.0
    while total_time<2.0:
        """ ODEによってVPythonのボールの位置を更新 """
        Ball_VS.pos =  Ball_Body.getPosition()

        vs.rate(1/dt)  # VPythonのアニメーションを進める
        world.step(dt)  # ODEの計算を進める
        total_time += dt
    vs.exit()

f:id:takeshiD:20160604125815g:plain

床をすり抜けてしまいました。
これは床とボールの衝突判定をしていないためすり抜けてしまいます。

コードを解説していきますが、その前にまずODEの概略について説明しましょう。

ODEの基本の2つ: 剛体と形状

ODEはまず、剛体(Body)と形状(Geom)の2つの要素から成ります。
剛体は質量特性を管理するオブジェクト
形状は衝突判定を管理するオブジェクトです

剛体も形状も、それぞれが一つのモノを表しますのでこれらを複数管理するオブジェクトが必要になります。
それがWorld(剛体の世界)とSpace(形状の世界)です。
以下に対応表を書いておきます。

種類	空間	含まれるオブジェクト
剛体・質量特性	World	Body
形状・衝突判定	Space	Geom

コードの解説~ODEの設定~

今回はODEの設定で

"""  ODE環境の設定  """

# ODEワールドの作成 + ワールドの重力設定
world = ode.World()
world.setGravity((0, -9.81, 0))

# 剛体の作成 + 質量特性の設定
Ball_Body = ode.Body(world) 
M = ode.Mass()              
M.setSphere(2500.0, 0.5)  
Ball_Body.setMass(M)        
Ball_Body.setPosition((0, 4, 0))

としてWorld(剛体の世界)しか作成していないので衝突判定ができなかったため、床をすり抜ける事態が起きました。
衝突判定についてはまた別の機会で上げます。

Worldの作成

剛体を計算させるためにはまず剛体が含まれる世界、Worldを作る必要があります。

# ODEワールドの作成 + ワールドの重力設定
world = ode.World()
world.setGravity((0, -9.81, 0))

ode.WorldでWorldのインスタンスworldを作成し、world.setGravityで重力を設定します。
この時、world.setGravityの引数は(x, y, z)という3次元のデータにしましょう。
タプル、リスト、VPythonのvector、numpyのndarrayのどれかにするといいでしょう。

VPythonではvector, ODEではタプルなどとするのは面倒なので、VPythonとODEの連携をするなら色々な高速な関数が入ってるnumpyを使うのがいいかもしれません。
今回は多くの例に習ってタプルにしました。

剛体の作成・配置

worldを作成したのでこの中に剛体を作成・配置していきましょう。

# 剛体の作成 + 質量特性の設定
Ball_Body = ode.Body(world) 
M = ode.Mass()              
M.setSphere(2500.0, 0.5)  
Ball_Body.setMass(M)        
Ball_Body.setPosition((0, 4, 0))

ode.Body(world)で、作成したworldにBody(剛体)を作成します。
剛体は質量特性を設定するオブジェクトなのでこのBodyにode.Massで質量特性を設定します。
M.setSphere(密度[kg/m³], 半径[m])としてMに球状の質量特性=密度+半径を設定します。
この場合、半径0.5[m]の質量分布は2500.0[kg/m³]で一定です。

次に作成した質量特性Mを剛体Ball_BodyにBall_Body.setMass(M)で設定します。
最後にBall_Body.setPosition*1で剛体の初期位置を設定します。

コードの解説~VPythonの設定~

VPythonの設定は完全に見た目だけのものです。
内部的な計算は全てODEでなされているのでVPythonではODEの設定をどう反映させるかだけが問題となります。

""" VPython環境の設定 """
Floor = vs.box(length=4, height=0.2, width=4)
Floor.pos = vs.vector(0, 0, 0)

Ball_VS = vs.sphere(pos=Ball_Body.getPosition(),
            radius=0.5,
                    color=vs.color.red)

床とボールの作成・設定

FloorとBall_VSを作成します。
vs.box(length, height, width)という形で引数を受け付けています。
(length, height, width)=(x, y, z)という形で軸が対応しています。

作成したFloorの位置をFloor.pos=vs.vector(x, y, z)という形で設定します。タプル、リスト、numpyのndarray、VPythonのvectorでも大丈夫です。今回はVPythonのvectorにしてみました。
この位置はboxの中心ではなく、6面のうちのどこかの面の中心が(x, y, z)に一致するようになっています。それだけではもちろんボックスの向きが分かりません。そこでFloor.axis=(x, y, z)という形でボックスの向きのベクトルを設定できます。

コードの解説~シミュレーションループ:VPythonとODEの連携~

実は連携というほど大げさなものじゃないです。

ODEで計算 → VPythonで見た目を描画

を繰り返すだけです。

if __name__=="__main__":
    dt = 0.05
    total_time = 0.0
    while total_time<2.0:
        """ ODEによってVPythonのボールの位置を更新 """
        Ball_VS.pos =  Ball_Body.getPosition()

        vs.rate(1/dt)  # VPythonのアニメーションを進める
        world.step(dt)  # ODEの計算を進める
        total_time += dt
    vs.exit()

dtはODEを計算するための時間ステップです。
dt=0.05[sec]ごとに計算が進むことになります。

Ball_VS.pos = Ball_Body.getPosition()によってODEで計算した結果をVPythonに渡します。

vs.rate(1/dt)によってVPythonの描画スピードを1/0.05=20[frames/sec]と設定します。
world.step(dt)によってODEの計算をdtだけ進めます。

whileループを外れたあと、vs.exit()によってVPythonの描画を終了しウィンドウを閉じます。

*1:0, 4, 0

2016-06-03

PyODEとVPythonのインストール

それぞれ特に難しくはありません。
私はWindowsなのでwhlファイルでインストールしました。 LinuxやMacの場合はpipで入れられるのかな？

PyODEのインストール
まず上記の公式サイト、ではなく！
Windows用Python拡張パッケージ非公式版(英語)
に飛んで、ページ内検索(CTRL+F)で「ode」と入力しましょう。
現在(2016/6/2)ではode-0.13.1が最新版です。
私はPython2.7 Windows7 64bitなのでode-0.13.1-cp27-none-win_amd64.whlをダウンロードしました。私はwhlファイルをよくインストールするのでC:\PythonPackageに一箇所に集めています。
今回も C:\PythonPackage\ode-0.13.1-cp27-none-win_amd64.whlとなるようにダウンロードしました。

[C:\PythonPackage]
$ ls
ode-0.13.1-cp27-none-win_amd64.whl

となっているので

[C:\PythonPackage]
$ pip install ode-0.13.1-cp27-none-win_amd64.whl

とすれば正常にインストールされます。
インストールされたか確認しておきましょう。

[C:\PythonPackage]
$ ipython
Python 2.7.10 (default, May 23 2015, 09:44:00) [MSC v.1500 64 bit (AMD64)]
Type "copyright", "credits" or "license" for more information.

IPython 4.1.1 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]: import ode

In [2]: world = ode.World()

worldを生成するときにエラーがでなければ大丈夫です。

VisualPythonのインストールまずVPythonWindows用ダウンロードページに飛びます。
あとはご自分のPythonが64bit, 32bitで適当な方を選んでバイナリをダウンロードしてダブルクリックすればインストーラーが勝手にインストールしてくれます。

2016-06-03

PyAudioの使い方まとめ

音声入出力ライブラリであるPyAudioの使い方をまとめた記事です。

PyAudioの基本メモ

PyAudioの基本メモ1 - たけし備忘録
 PyAudioの基本メモ2 音声入出力 - たけし備忘録

PyAudioとpyatgraphでリアルタイムプロット

pyqtgraph+PyAudioによるリアルタイムで音声プロット - たけし備忘録
pyqtgraph+PyAudioによるリアルタイムで音声プロットその2 - たけし備忘録

2016-06-02

PyODE+VisualPythonで物理シミュレーション目次

python VPython ODE 物理シミュレーション

一度はやってみたい物理シミュレーション
物理シミュレーションといえば物理演算エンジン色々な種類があります。
今回はその中でもPythonでバインディングされているC++の物理エンジンライブラリであるODE(OpenDynamicsEngine)を用いてシミュレーションをしてみます。

参考文献

ODE関連のサイト

DEMURA.NET(日本語)
C++ライブラリでのODEの丁寧な解説
ODEプログラミング解析so-zou.jp(日本語)
C++ライブラリでのODEの図ありの丁寧な解説
PyODE(英語)
ODEのPython バインディングであるPyODEの公式サイト
API Referenceをよく見る

VisualPython関連のサイト

VPython(英語)
VisualPythonの公式サイト
VisualモジュールVPython(日本語)
VPythonで主に使用する関数やオブジェクトの日本語解説　ありがたい
画像処理についてあれこれ: VPython
シンプルなコードが多く、辞書的でわかりやすい。テクスチャや光源の扱い方も載ってるのでとても助かる。

基本を学ぶ

1. PyODEとVPythonのインストール - たけし備忘録

2. ボールの自由落下 - たけし備忘録

3. ボールの自由落下(衝突判定有り) - たけし備忘録

2016-01-22

Julia+SymPyの勝手なまとめ

Julia python

自分の気になったところだけメモしていく形です
基本的にはJuliaのSymPyパッケージ内に入っているtutorial.mdというファイルから抜粋する形です。

SymPy の導入
シンボルの生成
- 1つずつシンボルを生成
- 複数のシンボルを一気に生成
特殊な定数 PI, E, oo
代入
シンボルから数値への変換
方程式の解
極限
- 基本の計算方法
- 右極限と左極限
数値計算

SymPy の導入

まずはjuliaを立ち上げてSymPyのパッケージをインストールしましょう。
PyCallに依存しているので入れてない方はPyCallもインストールしましょう。

$ julia
              _
   _       _ _(_)_     |  A fresh approach to technical co
  (_)     | (_) (_)    |  Documentation: http://docs.julia
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _  |  |
  | | |_| | | | (_| |  |  Version 0.4.3 (2016-01-12 21:37)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ r
|__/                   |  x86_64-w64-mingw32

julia> Pkg.add("PyCall")
...........なんかいろいろ出る
julia> Pkg.add("SymPy")
...........なんかいろいろ出る
julia>

<*追記*>
依存がPyCallとSymPyだけじゃありませんでした。
JuliaのPyCallとSymPyパッケージを入れた環境に、Python(2でも3でもいい)をインストールし、かつPythonのSymPyライブラリをインストールする必要があります。

python.jp
こちらのPython日本語ホームページからダウンロードできます。バージョンは2でも3でもお好きな方で大丈夫です。

Pythonのインストールが終わったら、次はPythonのSymPyをインストールする必要があります。
Pythonは2.7.9以上のバージョン(3を含む)ならば、PIPというパッケージ管理ツールが同梱していますのでそのままお使いいただけます。

[Linux, OSX, Windows]の場合

$ python -m pip install sympy

で基本的にはインストールできます。

上記でインストール出来ない場合(Windowsのみ記載)

上記の方法でインストール出来ない場合、次の非公式パッケージがダウンロードできるサイトで.whlというファイルをダウンロードしましょう。
Python非公式パッケージ
現時点(2016/1)ではsympy-0.7.6.1-py2.py3-none-any.whlが最新のようです。
この.whlファイルを、ターミナルから操作しやすいディレクトリに入れておきましょう。
今回の場合、C:\test\sympy-0.7.6.1-py2.py3-none-any.whlとして保存しておくことにします。

$ pwd # 現在地
C:\test

$ ls
sympy-0.7.6.1-py2.py3-none-any.whl

$ python -m pip install sympy-0.7.6.1-py2.py3-none-any.whl
(色々出る)

で完了です。

上記でインストール出来ないかつ不安だからpython2と3の両方インストールしちゃったよの場合(Windowsのみ記載)

pythonはデフォルトで2と3が共存できます(3の方が新しいですが2への後方互換性が無いです)。
そのような場合 pyコマンドで2と3を明示的に切り替えることが出来ます。

# Python2を使う場合
$ py -2 -m pip install sympy-0.7.6.1-py2.py3-none-any.whl

# Python3を使う場合
$ py -3 -m pip install sympy-0.7.6.1-py2.py3-none-any.whl

シンボルの生成

SymPyのJuliaでのラッパは色々なシンボルの生成の仕方があります。

1つずつシンボルを生成

julia> x = Sym("x") # Sym関数での生成
x
julia> @syms x      # symsマクロでの生成
(x,)
julia> x
x

複数のシンボルを一気に生成

julia> x, y, z = symbols("x, y, z")   # symbols関数での生成
(x,y,z)
julia> a, b, c, d, e = symbols("a:e") # ":"というイテレータで一気に生成
(a,b,c,d,e)
julia> @syms x y z                  # symsマクロで生成
(x,y,z)

マクロでの生成が一番手軽でよさそう

特殊な定数 PI, E, oo

PIは円周率
Eは自然定数の底(ネイピア数)
ooは無限大です

julia> PI,E,oo
(pi,E,oo)

juliaはデフォルトでpiという定数がありますがPIとは異なります。

julia> [asin(1), asin(Sym(1))]
2-element Array{SymPy.Sym,1}
[1.5707963267949]
[               ]
[      pi       ]
[      --       ]
[      2        ]

代入

シンボルで定義した式に、数値またはシンボルを代入します。
Julia固有の色々な代入の仕方があります。

subs関数による代入

前準備

julia> @syms x y
(x,y)

julia> f = x^2+2x+1
 2
x  + 2*x + 1

シンボルxにシンボルyを代入

julia> subs(f, x, y)
 2
y  + 2*y + 1

シンボルxに数値1を代入

julia> subs(f, x, 1)
4

便利な代入

前準備

julia> @syms x y z
(x,y,z)

julia> f=x+y+z
x + y + z

便利な代入1

これは便利な気がする

julia> f(x=>1, y=>pi)
z + 1 + pi

便利な代入2

代入の順序がよくわかんない

julia> f(1, pi)
y + 1 + pi

便利な代入3

パイプライン(|>という演算子)による代入です

julia> f |> subs(x, 1)
y + z + 1

便利な代入1の=>による代入が一番わかりやすい気がする。

シンボルから数値への変換

evalf関数とN関数の二つの方法があります。
NはBigFloatという型らしいです。よく知りません。
evalfはお馴染みです。

julia> evalf(PI)
3.14159265358979

julia> N(PI)
3.141592653589793

julia> evalf(PI, 30)            # 桁数指定 30桁
3.14159265358979323846264338328

julia> N(PI, 30)                # 桁数じゃなかった。何者だろうこの引数は
3.1415926535897932384626433832793

方程式の解

solve関数を使います。

前準備

julia> @syms x
(x,)

julia> f=x^2+x+1
 2
x  + x + 1

解を得る1 x^2+x+1=0

julia> solve(f)
2***element Array{SymPy.Sym,1}
[        ___  ]
[  1   \/ 3 *I]
[- - - -------]
[  2      2   ]
[             ]
[        ___  ]
[  1   \/ 3 *I]
[- - + -------]
[  2      2   ]

解を得る2 x^2+x+1=1

julia> solve(Eq(f, 1))
2-element Array{SymPy.Sym,1}
[-1]
[  ]
[0 ]

解を得る3 x^2+x+1=0

julia> solve(f⩵0)    
2-element Array{SymPy.Sym,1}
[        ___  ]
[  1   \/ 3 *I]
[- - - -------]
[  2      2   ]
[             ]
[        ___  ]
[  1   \/ 3 *I]
[- - - -------]
[  2      2   ]

# ※注意！！ f==0と書いていますがこれは \Equalの後にTabを押して出る特殊記号です。

極限

極限はlimit関数で操作できます。

基本の計算方法

julia> @syms x h
(x,h)

julia> limit(sin(x)/x, x, 0)
1

julia> limit(sin(x)/x, x=>0)
1

julia> limit((1+1/x)^x, x=>oo)
E

julia> limit((sin(x+h)-sin(x))/h, h=>0)
cos(x)

右極限と左極限

sign関数を使って例を書きます。
sign関数とは

$sign(x)=\left\{\begin{matrix} 1 & (x> 0) \\ 0 & (x=0) \\ ***1 & (x< 0 ) \end{matrix}\right.$

という関数です。つまり左極限と右極限が異なる不連続関数です。

julia> limit(sign(x), x=>0, dir="-") # 左極限
-1

julia> limit(sign(x), x=>0, dir="+") # 右極限
1

数値計算

やっぱりなんだかんだ数値計算が出来てなんぼさ。
ということで数値計算の主な方法です。

関数の定義

Juliaでの関数の定義はいろいろなやり方があります。
まずSymPyを使うとシンボルを定義します。
シンボルはSymPy.Sym型の変数ですのでJuliaの性質である多重ディスパッチが使えます。
多重ディスパッチとは、同名の関数、例えばtestという関数を

function test(x::Int)
    return x
end

function test(x::Float64)
    return 2x
end

として同名で、引数の型の違いによって出力や内容を変える性質のことです。

多重ディスパッチを使うと、

@syms x
f(x) = x^2+x+1

という風に定義すると

julia> f(x)
2
x  + 2*x + 1

julia> f(1)
4

julia> f(1.0)
4.0

julia> f(1:10)
ERROR: MethodError: `*` has no method matching *(::UnitRange{Int64}, ::UnitRange{Int64})
Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
  *{T}(::Bidiagonal{T}, ::AbstractArray{T,1})
  *(::Number, ::AbstractArray{T,N})
  ...
 in power_by_squaring at intfuncs.jl:78
 in f at none:1

1:10を引数にするとエラーが出ました。
そこでArray型にも対応できるように

f(x) = x.^2+2x+1

とします。 x^2がx.^2になっただけです。
このようにすれば

julia> f(1:10)
10-element Array{Int64,1}:
   4
   9
  16
  25
  36
  49
  64
  81
 100
 121

となります。
行列も引数に与えられます。

SymPyのままでも

@syms x
f = x.^2+2x+1

としても

julia> f(x)
 2
x  + 2*x + 1

julia> f(1)
4

julia> f(1.0)
4.00000000000000

julia> f(1:10)
ERROR: MethodError: `convert` has no method matching convert(::Type{SymPy.Sym}, ::UnitRange{Int64})
This may have arisen from a call to the constructor SymPy.Sym(...),
since type constructors fall back to convert methods.
Closest candidates are:
  SymPy.Sym(::Any)
  SymPy.Sym(::Any...)
  call{T}(::Type{T}, ::Any)
  ...
 in subs at C:\Users\takeshi\.julia\v0.4\SymPy\src\subs.jl:46
 in call at C:\Users\takeshi\.julia\v0.4\SymPy\src\subs.jl:115

とこのように関数のように使うことはできます。
しかし1:10などの配列や行列などは代入できません。
この理由としては
f = x.^2+2x+1
と定義するとSymPy.Sym型であり
f(x) = x.^2+2x+1
と定義するとこれは多重ディスパッチによって無名関数になるからです。
型を確認してみましょう。

区別がつきやすいように

@syms x
f(x) = x.^2+2x+1
g = x.^2+2x+1

としておきます。

julia> typeof(f)
Function

julia> typeof(g)
SymPy.Sym

fがFunction型、gがSymPy.Sym型になっています。

またさらにfの場合

julia> typeof(f(x))
SymPy.Sym

julia> typeof(f(1))
Int64

julia> typeof(f(1.0))
Float64

julia> typeof(f(1:10))
Array{Int64,1}

となります。
gの場合は

julia> typeof(g(x))
SymPy.Sym

julia> typeof(g(1))
SymPy.Sym

julia> typeof(g(1.0))
SymPy.Sym

となっていずれもSympy.Sym型になります。
なかなかこれが使えます。

微分

当たり前ですがPythonのSymPyと同様に記号微分が出来ます。

julia> @syms x
(x,)

julia> f(x) = x.^2+2x+1
f (generic function with 1 method)

julia> diff(f(x))
2*x + 2

julia> diff(f(x), x)
2*x + 2

julia> diff(f)
2*x + 2

diff関数は多重ディスパッチで
diff(f:SymPy.Sym)だけでなく
diff(f::Function)も定義されているので、無名関数として定義したfを入れても微分してくれます。

diff関数のオプション

diff関数は任意の回数だけ微分することができます。
以下に例を書きます。

julia> f(x)
 2
x  + 2*x + 1

julia> diff(f(x), x, 0)  # 0階微分
 2
x  + 2*x + 1

julia> diff(f(x), x, 1)  # 1階微分
2*x + 2

julia> diff(f(x), x, 2)  # 2階微分
2

もちろんfのままでも渡せます。が、少し文法が異なります。
要するにxが消えます。無名関数のまま渡すので、引数の情報がいらないわけです。

julia> f(x)
 2
x  + 2*x + 1

julia> diff(f, 0)  # 0階微分
 2
x  + 2*x + 1

julia> diff(f, 1)  # 1階微分
2*x + 2

julia> diff(f, 2)  # 2階微分
2

偏微分

2つのシンボルで定義された関数は偏微分できます。
さっそく例を出していきましょう。

julia> @syms x y
(x,y)

julia> f(x, y) = x^2 + y^2 + x*y
f (generic function with 2 methods)

julia> diff(f(x,y), x)    # xについて偏微分
2*x + y

julia> diff(f(x,y), y)    # yについて偏微分
x + 2*y

julia> diff(f(x,y), x, y) # x,yについて偏微分
1

julia> diff(f(x,y))
ERROR: PyError (:PyObject_Call) <type 'exceptions.ValueError'>
ValueError('\nSince there is more than one variable in the expression, the\nvariable(s) of differentiation must be supplied to differentiate x**2\n+ x*y + y**2',)
  File "C:\Python27\lib\site-packages\sympy\core\function.py", line 1638, in diff
    return Derivative(f, *symbols, **kwargs)
  File "C:\Python27\lib\site-packages\sympy\core\function.py", line 987, in __new__
    must be supplied to differentiate %s''' % expr))

 [inlined code] from C:\Users\takeshi\.julia\v0.4\PyCall\src\exception.jl:81
 in pycall at C:\Users\takeshi\.julia\v0.4\PyCall\src\PyCall.jl:361
 in fn at C:\Users\takeshi\.julia\v0.4\PyCall\src\conversions.jl:194
 in call_sympy_fun at C:\Users\takeshi\.julia\v0.4\SymPy\src\SymPy.jl:262
 in sympy_meth at C:\Users\takeshi\.julia\v0.4\SymPy\src\SymPy.jl:268
 in diff at C:\Users\takeshi\.julia\v0.4\SymPy\src\math.jl:190

julia> diff(f)    # 無名関数を偏微分
2*x + 2

積分

微分と同様に積分です。

不定積分

julia> f(x)
 2
x  + 2*x + 1

julia> integrate(f(x), x)
 3
x     2
-- + x  + x
3

julia> integrate(f(x))
 3
x     2
-- + x  + x
3

julia> integrate(f)
 3
x     2
-- + x  + x
3

より一般的な書き方

julia> @syms x n
(x,n)

julia> integrate(x^n, x)
/log(x)  for n = -1
|
| n + 1
<x
|------  otherwise
|n + 1
\

めちゃめちゃわかりにくいですがこれは

$\int x^n dx=\left\{\begin{matrix} \log(x) & (n=***1) \\ \frac{x^{n+1}}{n+1} & (n \neq***1) \end{matrix}\right.$

ということです。

定積分

julia> integrate(f, 0, 3)
21

julia> integrate(f, 0, x)
 3
x     2
-- + x  + x
3

julia> integrate(f, x, 3)
   3
  x     2
- -- - x  - x + 21
  3

julia> integrate(f, x, 2x)
   3
7*x       2
---- + 3*x  + x
 3

積分範囲もシンボルが使えます。

重積分

同様に重積分も行えます。

julia> @syms x y
(x,y)

julia> integrate(x*y, (y, 0, 1), (x, 0, 1))
1/4

julia> integrate(x^2*y, (y, 0, 1-x^2), (x, -1, 1))
8/105

シンボルを含んだ積分をした後に数値のみの積分も出来ます。

テイラー展開

テイラー展開もできちゃいます。
めんどくさいので略。

第n部分和と無限級数

基本的な使い方

まずは基本からです。summation関数を使います。

julia> @syms k n
(k,n)

julia> summation(2^k, (k, 1, n))    # 第n部分和
 n + 1
2      *** 2

julia> summation(1/k^2, (k, 1, oo)) # 無限級数(この例はゼータ関数)
  2
pi
---
 6

limit関数を使った無限級数

limitを使えば部分和を無限級数にできちゃいます。

julia> Sn = summation(1/k^2, (k, 1, n))
harmonic(n, 2)

julia> limit(Sn, n=>oo)
  2
pi
---
 6

harmonic(n, 2)はnをシンボルとした2次の調和級数です。要するにζ(2)です。