確率微分方程式のカオス展開による近似解法 -

適当な参考文献。

The asymptotic error of chaos expansion approximations for stochastic differential equations
https://arxiv.org/abs/1906.01209

Dynamical Polynomial Chaos Expansions and Long Time Evolution of Differential Equations with Random Forcing
https://arxiv.org/abs/1505.00047

Wiener Chaos Expansion and Numerical Solutions of Stochastic Partial Differential Equations
https://thesis.library.caltech.edu/1861/

多項式カオス展開と呼ばれる方法によって、確率常微分方程式を等価な無限個の常微分方程式に置き換えることが出来る。同様に、確率偏微分方程式も無限個の偏微分方程式に置き換えることが出来る。これらの置き換えは厳密に正しいが、無限個の微分方程式から、適当に有限個の方程式を選んで、近似計算に使うことができる。

確率微分方程式の数値計算への応用が現れたのは、1980年代か1990年代初頭あたりっぽい(ちゃんと調べてない)けど、多項式カオス自体は、Wienerの1938年の論文で導入された。
The Homogeneous Chaos
http://users.dimi.uniud.it/~giacomo.dellariccia/Table%20of%20contents/Wiener1938.pdf

論文タイトルは、"homogeneous chaos"だけど、本文中で、polynomial chaosと呼んでいる。現在、"カオス"という名前は、カオス理論で使われる方が有名になった。カオス理論と多項式カオス展開は、関係ない。ただ、Wienerは、上論文中で、

Physical theories of chaos, such as that of turbulence, or of the statistical theory of a gas or a liquid, may or may not be theories of equilibrium.

と書いているので、動機自体は、カオス理論と近い所にあったと思われる。個人的に、カオス理論は別に好きじゃないけど、カオス展開は面白いので、知名度が上がってほしい。

Wienerの論文は長い上に、少なくとも現代の人にとっては、あんまり読みやすくない。CameronとMartinの1947年の論文は短いし、現代の人でも読みやすい。

The Orthogonal Development of Non-Linear Functionals in Series of Fourier-Hermite Functionals
https://doi.org/10.2307/1969178

CameronとMartinの1947年の論文では、いわゆるCameron-Martin空間の正規直交基底がimplicitに導入されている。

標準Wiener過程の適当な説明。Wiener過程は、ホワイトノイズの積分として「定義」されるガウス過程であるという性質と、推移確率密度が熱核で書けるマルコフ過程であるという性質を持つ。2つの性質は、見かけ上異なるけど、どっちを出発点としてもいい。WienerがWiener過程を考えた時、まだマルコフ過程もガウス過程も定義されてなかったけど、Wienerの出発点は、どっちかというと後者に基づいていたっぽい(論文を見てないので、二次情報からの推測)。今となっては、前者の方が、Wiener過程が何故基本的なのか理解しやすいと思う。

確率解析でWiener過程が重要な理由も前者に由来すると思う。物理の教科書では、ホワイトノイズ $Z_t$ を、 $E[Z_t$ =0]かつ $E[Z_s Z_t]=\delta(s-t)$ を満たす"定常ガウス過程"もどきと見なして扱うのが一般的に行われている。数学では、自己相関がデルタ関数であるために、ホワイトノイズを直接扱いづらいので、ホワイトノイズの積分であるWiener過程をベースに理論を作るのが一般的な処方となった。

$W_t$ を標準Wiener過程とすると、形式的に
$W_{t} = \displaystyle \int_{0}^{t} Z_{s} ds$
と考えられる。Wiener過程の標本経路が微分不可能であることを思うと、ホワイトノイズのサンプルは、少なくとも、普通の関数であると考えることはできない。この積分を、正当化する方法は分からないけど、無視して進む。

$E[\cdot]$ で確率変数の期待値を表す時、期待値を取る操作が、和を取る操作と可換であることを、積分にも適用できると考えると
$E[W_{t}] = E \left[ \displaystyle \int_{0}^{t} Z_{s} ds \right] = \displaystyle \int_{0}^{t} E[Z_s] ds = 0$
と
$E[W_sW_t] = E \left[\displaystyle \int_{0}^{s} \int_{0}^{t} Z_u Z_v du dv \right] = \displaystyle \int_{0}^{s} \int_{0}^{t} E[Z_uZ_v] du dv = \displaystyle \int_{0}^{s} \int_{0}^{t} \delta(u-v) du dv = \mathrm{min}(s,t)$
が成立してほしい。この議論を正当化することはできない(か、出来ても却って準備が面倒になる)ので、この性質で特徴付けられる(非定常)ガウス過程というのが、標準Wiener過程の一つの定義になる。

$E[W_s W_t] = \mathrm{min}(s,t)$ が成立するなら、 $E[(W_t-W_s)^2]= t + s -2 \mathrm{min}(s,t)=|t-s|$ などから、増分が分散 $|t-s|$ の正規分布に従うと期待される。また、 $t_1 \lt t_2 \lt t_3$ に対して、 $W_{t_3}-W_{t_2}$ と $W_{t_2}-W_{t_1}$ は、
$E[(W_{t_3} - W_{t_2})(W_{t_2}-W_{t_1})] = 0$
で互いに直交し、正規分布の性質から、重なりのない区間の増分は互いに独立であることも分かる。推移確率密度が熱核で書けるマルコフ過程になっているという出発点からも、同じ要請が出る。

$h \in L^2[0,1]$ に対して、 $\displaystyle \int_{0}^{1} h(t) dW_t$ が、平均0で分散が $h$ のノルムの二乗で与えられる正規分布に従うことは、積分の定義に戻れば、すぐに理解できる( $h$ は確定的な関数なので、特に難しいことはない)

$L^2[0,1]$ の正規直交基底 $e_n(n=1,2,\cdots)$ を適当に固定して、 $\chi_I(t)$ は、 $t \in I$ の時に1、それ以外では0となる関数とする。 $0 \le s \le 1$ に対して、
$c_n = \displaystyle \int_{0}^{1} \chi_{[0,s]}(t) e_n(t) = \displaystyle \int_{0}^{s} e_n(t)$
を定義して、
$W_s = \displaystyle \int_{0}^{s} dW_{t} = \displaystyle \int_{0}^{1} \chi_{[0,s]}(t)dW_{t} = \displaystyle \int_{0}^{1} \sum_{n=1}^{\infty} c_n e_n(t) dW_{t}$
で、形式的には
$W_s = \displaystyle \sum_{n=1}^{\infty} \left(\displaystyle \int_{0}^{1} e_n(t) dW_{t} \right) c_n = \displaystyle \sum_{n=1}^{\infty} \left(\displaystyle \int_{0}^{1} e_n(t) d W_{t} \right) \displaystyle \int_{0}^{s} e_n(t)$
という式変形ができる。実用上は、どうせ有限和しか取らないけど、収束するかどうかについては、もう少し議論がいる。

証明は、例えば、1966年のSheppという人の論文に書かれていて(Theorem4)、”This result appears to be new except for two special cases of Wiener and Levy"と述べられている。
Radon-Nikodym derivatives of gaussian measures
https://doi.org/10.1214/aoms/1177699516

どうやって証明するのか知らないけど、Kolmogorovのthree-series theoremというものがあるらしく、それを使っている。Kolmmogorovが1925年に証明して、確率論の人には、よく知られた定理らしい。原理的には、勝手な正規直交基底を取っていいとしても、数値計算する時には、有限項で打ち切ることになるので、誤差評価ができることが望ましい。

一般のGauss過程に対しては、類似の定理として、Karhunen–Loèveの定理と呼ばれるものがある(Wiener過程に対しては、基底の形が制限されるので特殊ケースしか示せないけど)。この定理は、1947年に、定常ガウス過程に対して、KacとSiegertによって述べられている。

An explicit representation of stationary gaussian process
https://doi.org/10.1214/aoms/1177730391

KacとSiegertが定常性を仮定した理由は、分からないけど、(多分ガウス過程を最初に導入したと思われる)1934年のKhintchineの論文でも、定常性が仮定されているので、それに従っただけかもしれない。

Korrelationstheorie der stationären stochastischen Prozesse
https://doi.org/10.1007/BF01449156

KarhunenとLoèveも、1947年前後、独立に、なにか論文を書いたらしいけど、電子化されたものを発見できなかったので内容の詳細は不明。1951年のCramérの論説に2人への言及が見られる。
A Contribution to the Theory of Stochastic Processes
https://projecteuclid.org/proceedings/berkeley-symposium-on-mathematical-statistics-and-probability/Proceedings-of-the-Second-Berkeley-Symposium-on-Mathematical-Statistics-and/Chapter/A-Contribution-to-the-Theory-of-Stochastic-Processes/bsmsp/1200500237

何を典拠にするかはともかく、平均0、分散1の正規分布に従う無限個の独立な確率変数 $\xi_i$ ( $i=1,2,\cdots$ )と $L^2[0,1]$ の正規直交基底 $e_n$ を適当に取って、標準Wiener過程を
$W_t= \displaystyle \sum_{n} \xi_n \displaystyle \int_{0}^{t} e_{n}(u)du$
と表示することができる。上の話は全部忘れて、これを定義のように思っても、差し支えない。

$\xi_1,\cdots,\xi_n$ は互いに独立で、標準正規分布に従うので、
$E[f(\xi_1 , \cdots , \xi_n)] = (2 \pi)^{-n/2} \displaystyle \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} f(u_1,\cdots ,u_n) e^{-(u_1^2 + \cdots + u_n^2)/2} du_1 \cdots du_n$
が成り立つのは、直感的に明らかで、 $\xi_1,\cdots,\xi_n$ の関数を、Hermite多項式の直積を基底にして展開したくなる。これが、polynomial chaosの"polynomial"の起源。

多重指数 $\alpha=(\alpha_1,\alpha_2,\cdots)$ ( $\alpha_j$ は自然数で $\sum_{n} \alpha_n \lt \infty$ )に対して
$T_{\alpha}(\mathbf{\xi}) = \displaystyle \prod_{n=1}^{\infty} H_{\alpha_n}(\xi_n)$
を定義する。 $H_k(x)$ は正規化したHermite多項式。

一般の標準Wiener過程の汎関数は
$f[W_t] = \displaystyle \sum_{\alpha} c_{\alpha} T_{\alpha}(\mathbf{\xi})$
という形で展開できる。この展開は、 $L^2[0,1]$ の正規直交基底 $e_n$ の取り方に依存していることは忘れるべきでない。

ここで、
$c = \displaystyle \sum_{\alpha} T_{\alpha}(\mathbf{\xi})$
に対して
$E[c T_{\beta}(\mathbf{\xi})] = c_{\beta}$
なので、 $c=0$ という確率論的な条件は、 $T_{\alpha}$ の係数が全部0という条件と同じ。

確率常微分方程式の場合は、係数 $c_{\alpha}$ が時間に依存し、確率偏微分方程式なら位置(や時間)に依存すると考える。係数自体は、普通の(確率的でない)微分方程式に従い、全ての係数を決定すれば、確率微分方程式が解けたことになる。

簡単な例として、Ornstein-Uhlenbeck過程の従う確率常微分方程式
$d x_{t} = (a - b x_t) dt + c \cdot dW_{t}$
$x_{t=0} = x_{0}$
を考える。これは、Langevin方程式を少し一般化したものと見なせる。 $x_t$ は、Langevin方程式の文脈では、時間tにおける粒子の速度と解釈される。物理の教科書では、ホワイトノイズ $Z_{t}$ が定義できるとして、 $dW_{t} = Z_{t} dt$ であるかのように書かれてるのが普通。

Langevinは1908年に、(フランス語タイトルを英訳すると)"On the Theory of the Brownian Motion"という論文を書き、UhlenbeckとOrnsteinは1930年に、"On the Theory of the Brownian Motion"という論文を書いている。
Sur la théorie du mouvement brownien
https://fr.m.wikisource.org/wiki/Sur_la_th%C3%A9orie_du_mouvement_brownien

On the Theory of the Brownian Motion
https://doi.org/10.1103/PhysRev.36.823

論文を読むと、方程式をLangevin方程式と呼び、その解をOrnstein-Uhlenbeck過程と呼ぶのは、一応、妥当に思える(他に、貢献者がいないかどうかは分からないが)。Doobは1942年の論文で、Ornstein-Uhlenbeck過程を、(非自明な)定常Gauss-Markov過程として特徴付けた。
The Brownian Movement and Stochastic Equations
https://doi.org/10.2307/1968873

重要というわけではないけど、以下の論文は、読み物として面白い。
Laplace and the origin of the Ornstein-Uhlenbeck process
https://doi.org/10.2307/3318524

で、変数 $x_t$ を
$x_t = \displaystyle \sum_{\alpha} x_{\alpha}(t) T_{\alpha}(\mathbf{\xi})$
と展開する。

形式的に計算して
$d x_{t} = \displaystyle \sum_{\alpha} x_{\alpha}'(t) T_{\alpha}(\mathbf{\xi}) dt$
$(a - b x_{t})dt = (a - b \displaystyle \sum_{\alpha} x_{\alpha}(t) T_{\alpha}(\mathbf{\xi})) dt$
$c \cdot dW_{t} = \displaystyle \sum_{\alpha} \xi_{n} e_{n}(t) dt$
と思うと、無限個の微分方程式
$x_{\alpha}'(t) = (a I_{\{\alpha=0\}} - b x_{\alpha}(t)) + c \displaystyle \sum_{n} e_{n}(t) I_{\{\alpha=\epsilon_n\}}$
$x_{\alpha}(0) = I_{\{\alpha=0\}} x_0$
に書き直すことが出来る。

$\epsilon_n$ は、k番目だけ1で、他は0という形の多重指数( $(\epsilon_n)_{k}=\delta_{nk}$ )。従って、 $T_{\epsilon_n}(\mathbf{\xi}) = H_1(\xi_n) = \xi_n$ となる。 $I_{\{\alpha=\beta\}}$ は、 $\alpha=\beta$ の時、1で、それ以外の時は0となる定数。

一般には、確率変数の積を普通の積と同じように扱って、 $F \cdot dW_{t} = (F \cdot Z_{t})) dt$ みたいな形式的計算をやると、間違った結果が得られる。Ornstein-Uhlenbeck過程の場合には、確率変数の積が出てこないので問題なかった。しかし、例えば、幾何ブラウン運動の満たす確率微分方程式だと、失敗する。それでも、確率変数の積をWick積で解釈してやると、うまくいくっぽい( $f(t, x_t)dt$ みたいな式を計算する時は、確率変数の積も、普通の積として計算していい)。何でWick積にするといいのか、直感的な説明が欲しいけど分からない。

cf)幾何ブラウン運動で、うまくいかない理由を検討した時に、
chaos expansion methods in malliavin calculus. a survey of recent results
https://sites.dmi.uns.ac.rs/nsjom/Papers/45_1/NSJOM_45_1_045_103.pdf
のRemark5.9とかを見た。他に、
ホワイトノイズ解析の新展開
https://www.math.is.tohoku.ac.jp/~obata/research/file/2008-IIAS-report.pdf
の5.3節とか。しかし、何で、Wick積にするのかは、(私にとっては)magicのままである。

得られた常微分方程式は、3パターンに場合分けできて、 $\alpha=\mathbf{0}$ の時は
$x_{\mathbf{0}}'(t) = a - b x_{\mathbf{0}}(t)$
$x_{\mathbf{0}}(0) = x_0$
となり、解は
$x_{\mathbf{0}}(t) = x_0 e^{-bt} + \dfrac{a}{b}(1 - e^{-bt})$

第2のパターンは、 $\alpha=\epsilon_n$ の時で
$x_{\epsilon_n}'(t) = -b x_{\epsilon_n}(t) + c \cdot e_n(t)$
$x_{\epsilon_n}(0)=0$
となり、解は
$x_{\epsilon_n}(t) = c \cdot e^{-bt} \displaystyle \int_{0}^{t} e_n(u) e^{bu} du$

第3のパターンは、 $|\alpha| = \sum_{n=1}^{\infty} \alpha_n \ge 2$ の時で
$x_{\alpha}'(t) = -b x_{\alpha}(t)$
$x_{\alpha}(0) = 0$
で、解は
$x_{\alpha}(t)=0$

一方、確率積分を使って解を書けば
$x_t = x_0 e^{-bt} + \dfrac{a}{b}(1-e^{-bt}) + c \cdot e^{-bt} \displaystyle \int_{0}^{t} e^{bs} \cdot dW_s$
で、両者は同じ結果を与える。

Ornstein-Uhlenbeck過程に対しては、全ての係数が満たす微分方程式を解けるので、カオス展開は、厳密解を与える。一般的には、カオス展開で与えられる微分方程式が解けるとは限らない。けど、係数の式の形が単純な場合は、微分方程式そのものは書くことができる。なので、(確率的でない)微分方程式を数値的に解くことで近似解が得られる。そこでは、伝統的な微分方程式の数値計算アルゴリズムが、自由に使える。

カオス展開で得られた微分方程式の解が分かれば、適当な個数の標準正規分布に従う乱数を振って、エルミート多項式に代入して、
$\displaystyle \sum x_{\alpha}(t) T_{\alpha} (\mathbf{\xi})$
を計算するだけで、一つのpath(の近似)をサンプリングできる。

確率常微分方程式の数値解法として、Euler-丸山法というのが、よく出てくる。これは、常微分方程式のEuler法を、確率常備分方程式に、そのまま拡張したもので、一番簡単なアルゴリズムと思う。
Euler–Maruyama method
https://en.wikipedia.org/wiki/Euler%E2%80%93Maruyama_method

Ornstein-Uhlenbeck過程を例にして、カオス展開による数値解と、Euler-丸山法による解と比較することを考える。カオス展開での具体的な計算には、 $L^2[0,1]$ の正規直交基底を決める必要があるので、
$e_{1}(t) = 1$
$e_{n}(t) = \sqrt{2} \cos( (n-1)\pi t)$ $(n=2,3,\cdots)$
と取る。これは、Wienerが使ったらしい基底。

Kを適当に選んだ自然数として、正規乱数 $\xi_1,\cdots,\xi_{K}$ に対して、近似Wiener過程を、以下の有限和で生成する。
$W(t) = \displaystyle \sum_{n=1}^{K} \xi_n \displaystyle \int_{0}^{t} e_n(u)du$

実験した限り、 $t \gt s$ に対して、 $W(t)-W(s)$ は、 $t-s$ がある程度大きい時には、分散 $t-s$ の正規分布に従ってそうな感じがする。 $t-s$ を、小さくすると、期待される分布からの乖離は大きくなって、分散の大きさが期待される値と10倍位違うということも起こる。

$t-s$ をどこまで小さく出来るかは、Kの大きさに依存しているように思われる。標準正規乱数 $X$ を使って、 $W(t)=W(S) + \sqrt{|t-s|}X$ によって逐次的にWiener過程を生成する場合でも、|t-s|より小さい時間幅では線形補間されるだけので、これが悪いわけではない。結局、近似の代償は、同じ部分に出るということらしい。

近似Wiener過程のサンプルを繰り返し作って、それを $W^{(n)}(t)$ として、
$E[W_s W_t] \approx \dfrac{1}{N} \displaystyle \sum_{n=1}^{N} W^{(n)}(s) W^{(n)}(t)$
を評価すると、 $min(s,t)$ に近い値が得られる。こうした実験から、カオス展開によるWiener過程の近似は、悪くなさそうだと分かる。

近似Wiener過程を使って、Euler-丸山法では、Ornstein-Uhlenbeck過程は
$Y(0) = x_0$
$Y(t_{n+1}) = Y(t_{n}) + (a - b \cdot Y(t_n))(t_{n+1} - t_n) + c(W(t_{n+1}) - W(t_{n})$
で計算される。近似Wiener過程のサンプルごとに、近似Ornstein-Uhlenbeck過程のサンプルが一つ作れる。

通常は、正規乱数 $X$ を使って、逐次的に $W(t_{n+1}) = W(t_{n}) + \sqrt{t_{n+1} - t_{n}} X$ とWiener過程のサンプリングを行っていく。カオス展開で作った近似Wiener過程を使う場合は、 $t_{n+1}-t_{n}$ が小さすぎると、 $W(t_{n+1}) = W(t_{n})$ の分布が期待したものから大きく外れてるので、時間幅を小さく取りすぎないように注意する必要があると思われる。

一方、カオス展開では、Ornstein-Uhelenbeck過程は、
$X(t) = x_0 e^{-bt} + \dfrac{a}{b}(1-e^{-bt}) + c \displaystyle \sum_{n=1}^{K} \xi_{n} \displaystyle \int_{0}^{t} e_n(u) e^{b(u-t)} du$
で近似される。近似Wiener過程のサンプルを作るのに使用したのと同じ乱数を使って、近似Ornstein-Uhlenbeck過程のサンプルを作ると、Euler-丸山法による近似解の比較が出来る。

積分の計算が怠いけど、数値計算には必要なのでサボらずにやれば以下の通り。
$X(t) = x_0 e^{-bt} + \dfrac{a}{b}(1-e^{-bt}) + c \xi_{1} \dfrac{1 - e^{-bt}}{b} + \sqrt{2} c \displaystyle \sum_{n=1}^{K-1} \xi_{n+1} \dfrac{(b \cos(\pi n t) + \pi n \sin(\pi n t)) - b e^{-bt}}{b^2+\pi^2n^2}$

カオス展開、dt=1/50のEuler-丸山法、dt=1/500のEuler-丸山法の3つで比較したのが以下のグラフ。タイムステップを細かくしたEuler-丸山法と、Kを大きくしたカオス展開の結果は、ほぼ一致している。まぁ、そもそも、Kを有限に取ったら、確率的でない微分方程式を考えるのと、何ら変わらないので、当たり前という気もするけど、とりあえず、積分の計算が間違ってなさそうというくらいは分かる。

f:id:m-a-o:20211114233000p:plain

import math
import random
import numpy as np
import matplotlib.pyplot as plt

class WienerPath:
    def __init__(self, _rnds):
        self.zs = _rnds
    def __call__(self, t):
        assert(t>=0 and t<=1)
        K = len(self.zs)
        fvals = [t]
        fvals.extend([math.sqrt(2)*math.sin(k*math.pi*t)/(k*math.pi) for k in range(1,K)])
        ret = 0
        for (z,c) in zip(self.zs , fvals):
           ret += z*c
        return ret


#-- Ornstein-Uhlenbeck process by chaos expansion
class XPath:
    def __init__(self , _rnds , a=1.05 , b=0.7 , c=0.5 , y0=0):
        self.zs = _rnds
        self.a = a
        self.b = b
        self.c = c
        self.y0 = y0
    def __call__(self , t):
        assert(t>=0 and t<=1)
        K = len(self.zs)
        a,b,c=self.a,self.b,self.c
        fvals =  [c*(1-math.exp(-b*t))/b]
        ret = (self.y0)*math.exp(-b*t) + a*(1-math.exp(-b*t))/b
        fvals.extend([math.sqrt(2)*c*(b*math.cos(math.pi*n*t) + math.pi*n*math.sin(math.pi*n*t) -b *math.exp(-b*t))/(b*b+(math.pi*n)**2) for n in range(1,K)])
        for (z,c) in zip(self.zs , fvals):
           ret += z*c
        return ret


#-- Ornstein-Uhlenbeck process by Euler-Maruyama method
class YPath:
    def __init__(self , _W , a=1.05 ,b = 0.7 , c=0.5 , y0=0):
        self.W = _W
        self.a = a
        self.b = b
        self.c = c
        self.y0 = y0
    def trace(self , ts):
        t_n = 0
        y = self.y0
        ys = []
        a,b,c=self.a,self.b,self.c
        for t in ts:
           if t <= t_n:
              ys.append( y )
           elif t>1:
              pass
           else:
              y_next = y + (a-b*y)*(t - t_n) + c*(W(t) - W(t_n))
              t_n = t
              ys.append(y_next)
              y = y_next
        return np.array(ys)



if __name__=="__main__":
   random.seed(1934)   #-- for experiment
   W = WienerPath([random.gauss(0,1) for _ in range(5000)])
   X = XPath(W.zs)
   X50 = XPath(W.zs[:50])
   X500 = XPath(W.zs[:500])
   Y = YPath(W)
   ts = np.linspace(0,1,num=51)
   ts2 = np.linspace(0,1,num=501)
   ts3 = np.linspace(0,1,num=2001)
   ys = Y.trace(ts)
   plt.plot(ts3,np.vectorize(X)(ts3), label="chaos expansion (K=5000)" , linewidth=3 , alpha=0.4, color="blue")
   plt.plot(ts3,np.vectorize(X500)(ts3), label="chaos expansion (K=500)" , linewidth=3, color="pink" , alpha=0.8, linestyle="dashed")
   plt.plot(ts3,np.vectorize(X50)(ts3), label="chaos expansion (K=50)" , linewidth=3, color="black" , linestyle="dashed")
   plt.plot(ts,ys , label="Euler-Maruyama(dt=1/50)" , linestyle="dashed" , linewidth=3 , color="red")
   plt.plot(ts2,Y.trace(ts2), label="Euler-Maruyama(dt=1/500)" , linestyle="dotted" , linewidth=4  , color="grey")
   plt.legend()
   plt.show()