2011-08-23

数値計算

計算

非線形連立方程式
解の数値計算の方法
こちらのpdf
- $f_i(\vec{X}+\Delta \vec{X}) = f_i(\vec{X}) +\sum \frac{\par f_i}{\par x_j} \Delta x_j + O(|\Delta \vec{X}| ^2)$ より
- $\vec{f}(\vec{X}+\Delta \vec{X}) = \vec{0$ となるように $\Delta{\vec{X}}$ を選ぶとすると、2次より高次の項を無視すると
- $J \Delta \vec{X} = - \vec{f}(\vec{X})$ を得る
- $\vec{X}_{new} = \vec{X}_{old}+\Delta\vec{X}$ を繰り返し計算する
以前の rootSolve の記事

f<-function(x,y){(x-3)^2+y^2-3}
g<-function(x,y){sin(x)+exp(y-1)-1}

# fx<-function(x){2*(x-3)}
# fy<-function(y){2*y}
# gx<-function(x){cos(x)}
# gy<-function(y){exp(y-1)}

x<-0
y<-0
T<-100
d<-0.01

for(t in 1:T){

A<-matrix(c(
(f(x+d,y)-f(x-d,y))/(2*d),(f(x,y+d)-f(x,y-d))/(2*d),
(g(x+d,y)-g(x-d,y))/(2*d),(g(x,y+d)-g(x,y-d))/(2*d)
),byrow=T,2,2)


dX<- -solve(A)%*%c(f(x,y),g(x,y))
x<-x+dX[1]
y<-y+dX[2]
}
print(x);print(y)

実行

# 実行
> print(x);print(y)
[1] 1.997356
[1] -1.41234

2011-08-16

Solving Equation on the Computer

シミュレーション計算

Solving Equation on the Computer とこちらのコメントで教えていただいた数値計算について
数値計算
- Euler's method
- Improved Euler method
- Runge-Kutta method
- その他（Wiki）
$\dot{x} = f(x)$ について
Euler's method
- $x_{n+1}=x_{n}+f(x_n)\Delta t$
Improved Euler method
- 2点の傾きの平均をとる
- $\tilde{x}_{n+1} = x_{n}+f(x_n)\Delta t$
- $x_{n+1}=x_n+\frac{f(x_n) + f(\tilde{x}_{n+1})}{2}\Delta t$
Runge-Kutta method
- $\begin{matrix} k_1 &=& f(x_n)\Delta t \\ k_2 &=& f(x_n+\frac{k_1}{2})\Delta t \\k_3 &=& f(x_n+\frac{k_2}{2})\Delta t\\ k_4 &=& f(x_n+k_3)\Delta t \end{matrix}$
- $x_{n+1} = x_n + \frac{k_1+2k_2+2k_3+k4}{6}$

ソースで用いた条件は以下
- $\dot{x}=e^{-x}$
- 初期条件： $x(0)=0$
- $0 \leq t \leq T$ の範囲の計算

f<-function(x){exp(-x)}
T0<-0
T<-10
x0<-0

# 解析解
F<-function(t){log(t+1)}
t<-seq(from=T0,to=T,by=0.01)
plot(t,F(t),type="l",col=1,ylim=c(0,F(T)))

# オイラー法
dt<-0.5
Nt<-T/dt
xold<-x0
xn<-c(x0)
for(nt in 1:Nt){
	xnew<-xold+f(xold)*dt
	xn<-c(xn,xnew)
	xold<-xnew
	}

t<-seq(from=T0,to=T,by=dt)
par(new=T)
plot(t,xn,col=2,ylim=c(0,F(T)))

# improved-オイラー法
dt<-0.5
Nt<-T/dt
xold<-x0
xn<-c(x0)
for(nt in 1:Nt){
	x1<-xold+f(xold)*dt
	xnew<-xold+(f(xold)+f(x1))/2*dt
	xn<-c(xn,xnew)
	xold<-xnew
	}

t<-seq(from=T0,to=T,by=dt)
par(new=T)
plot(t,xn,col=3,ylim=c(0,F(T)))

# ルンゲ-クッタ法：fourth-order Runge-Kutta
dt<-1
Nt<-T/dt
xold<-x0
xn<-c(x0)
for(nt in 1:Nt){
	k1<-f(xold)*dt
	k2<-f(xold+k1/2)*dt
	k3<-f(xold+k2/2)*dt
	k4<-f(xold+k3)*dt
	xnew<-xold+(k1+k2*2+k3*2+k4)/6
	xn<-c(xn,xnew)
	xold<-xnew
	}

t<-seq(from=T0,to=T,by=dt)
par(new=T)
plot(t,xn,col=4,ylim=c(0,F(T)))

legend(x=c(4,10),y=c(0.5,1),c("Analytical solution","Euler's method","improved Euler method","fourth-order Runge-Kutta method"),col=c(1,2,3,4),lty=c(1,0,0,0),pch=c(-1,1,1,1))

2011-08-15

場の成分の計算

シミュレーション

ポテンシャルから場を計算する
- こちらの2.7 にもポテンシャルがでてくる
つまりは一階偏微分係数を求める

Nx<-10
Ny<-10

V<-matrix((1:(Nx*Ny))^2,Nx,Ny,byrow=T)
V_x<-matrix(0,Nx,Ny)

V_x[c(-1,-Nx),]<-(V[c(-1,-2),]-V[c(-(Nx-1),-Nx),])/2
V_x[1,]<-V[2,]-V[1,]
V_x[Nx,]<-V[Nx,]-V[Nx-1,]

persp(V,col=4)
persp(V_x,col=4)

# gx<-gxy<-1
# Ix<-matrix(0,Nx-1,Ny)
# Ix<-(-1)*gx*(V[-1,]-V[-Nx,])+(-1)*gxy*(V_x[-1,]+V[-Nx,])/2

2011-08-15

拡散

3次元の場合
は次のように表されている
- $\begin{cases} x &=& r \sin \phi \cos \theta \\ y &=& r \sin \phi \sin \theta \\ z &=& r \cos \phi \end{cases}$
については以下の関係も得られる
- $\tan \theta = \frac{y}{ x} ,\,\tan \phi = \frac{\sqrt{ x^2 + y^2}}{ z}$
このときを次のように変形できる
- $\begin{matrix} \frac{\par^2}{\par r^2} &=& (\frac{\par}{\par x} \frac{\par x}{\par r} + \frac{\par}{\par y} \frac{\par y}{\par r} +\frac{\par}{\par z} \frac{\par z}{\par r})^2 \\ &=& \frac{\par^2 }{\par x^2} (\frac{\par x}{\par r})^2 + \frac{\par^2}{\par y^2} (\frac{\par y}{\par r})^2 +\frac{\par^2}{\par z^2}(\frac{\par z}{\par r})^2 +(\frac{\par^2}{\par x \par y}+\frac{\par^2}{\par y \par x})\frac{\par x}{\par r}\frac{\par y}{\par r} + (\frac{\par^2}{\par y \par z}+\frac{\par^2}{\par z \par y})\frac{\par y}{\par r}\frac{\par z}{\par r}+(\frac{\par^2}{\par z \par x}+\frac{\par^2}{\par x }\par z)\frac{\par z}{\par r}\frac{\par x}{\par r} \end{matrix}$
たとえばのときのを計算してみると
- $\sin \phi \cos \phi \sin \theta = \frac{1}{2}$
- などとなっている

2011-08-14

Flows on the Line

力学

2. Flows on the Line
- 2.0 Introduction
- 2.1 A Geometric Way of Thinking
- 2.2 Fixed Points and Stability
- 2.3 Population Grouth
- 2.4 Linear Stability Analysis
- 2.5 Existence and Uniqueness
- 2.6 Impossibility of Oscillations
- 2.7 Potentials
- 2.8 Solving Equation on the Computer

2.0 Introduction
- $\dot{x} = f(x)$ の一次元(one-dimentional, first-order system)の挙動を扱う
- ここで関数 $f(x)$ は $x(t)$ の関数としておく
- $f(x,t)$ ではsecond-order になるのでここでは扱わない
2.1 A Geometric Way of Thinking
- $f(x)$ のグラフを $x-\dot{x}$ 平面に書くことでベクトル場を描き、安定解、不安定解などを判定する
2.2 Fixed Points and Stability
- グラフを描いて、平衡解とその安定性を判定する
2.3 Population Growth
- $f(x)$ が与えられ、そのグラフから解だけでなく時間発展の様子を判断する
2.4 Linear Stability Analysis
- 平衡解 $x^*$ の近くでの挙動を線型の方程式で近似する
- 小さい揺らぎ $\eta(t) = x(t) - x^*$ を与える
- $\dot{\eta} = \dot{x} = f(x) = f(\eta + x^*) = \sum_{n=0}^\infty \frac{f(x^*)}{n!} \eta^n = f(x^*) +f^{(1)}(x^*) \eta + O(\eta^2)$ により
- で与えられる
  - $f'(x*) > 0$ で揺らぎが増加、 $f'(x*) < 0$ で揺らぎが減衰する
  - 時定数も計算できる
  - $f'(x^*) =0$ の時は、Tayler展開の高次の項が関わる
2.5 Existence and Uniqueness
- 初期値が不安定解の場合、解が一意でないことがある
- Existens and Uniquness Theorem
  - 初期値問題について実数軸の開区間上でが連続でが実軸上の点ならば、初期値問題はある区間に解が存在し、一意である。
    - つまり、 $f(x)$ が十分滑らかならば解が存在して一意になるが、任意の時間で解が存在するという訳ではない
2.6 Impossibility of Oscillations
- これまでの $\dot{x} = f(x)$ の例では振動がみられない
- 非常に粘性の高い溶液中のバネによる運動のモデルもこれに相当する $m\ddot{x} \ll b\dot{x}$
2.7 Potentials
- ポテンシャルの導入による安定性の判定
- $f(x) = - \frac{dV}{dx}$ によりポテンシャルを導入する
- ここで $\frac{dV}{dt} = \frac{dV}{dx}\frac{dx}{dt} =- \{ \frac{dx}{dt} \} ^2 \leq 0$ なのでポテンシャルは時間とともに減少する
- ポテンシャルの極値をあたえるが平衡解で、それが下に凸（極小）ならば安定である
  - 平衡解について $\frac{dV}{dx} = 0$ と $f(x) = 0$ が同値
  - 安定性について $\frac{d^2V}{dx^2} > 0$ と $f' (x) <0$ が同値なので
2.8 Solving Equation on the Computer

2011-08-14

力学の教科書

力学

Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Studies in Nonlinearity)

作者: Steven H. Strogatz
出版社/メーカー: CRC Press
発売日: 2001/01/15
メディア: ペーパーバック
購入: 2人クリック: 46回
この商品を含むブログ (5件) を見る

- 1. Overview
Part 1. One-Dimensional Flows
- 2. Flows on the Line
  - 2.0 Introduction
  - 2.1 A Geometric Way of Thinking
  - 2.2 Fixed Points and Stability
  - 2.3 Population Grouth
  - 2.4 Linear Stability Analysis
  - 2.5 Existence and Uniqueness
  - 2.6 Impossibility of Oscillations
  - 2.7 Potentials
  - 2.8 Solving Equation on the Computer
- 3. Bifurcations
- 4. Flows on the Circle
Part 2. Two-Dimensional Flows
- 5. Linear Systems
- 6. Phase Plane
- 7. Limit Cycles
- 8. Bifurcations Revisited
Part 3. CHAOS
- 9. Lorenz Equations
- 10. One-Dimensional Maps
- 11. Fractals
- 12. Strange Attractors

2011-08-12

拡散

メモ
3次元はこちらへ
$\begin{matrix}\frac{\par^2}{\par r^2} & = & \frac{\par^2}{\par x^2} \cos^2\theta + \frac{\par^2}{\par y \par x}\sin \theta \cos \theta +\frac{\par^2}{\par x \par y} \sin \theta \cos \theta + \frac{\par^2}{\par y^2} \sin^2 \theta \\ & = & ( \frac{\par}{\par x} \cos\theta +\frac{\par}{\par y} \sin\theta)^2 \end{matrix}$
について
- $\theta = \frac{\pi}{4}$ のとき $\frac{1}{2}$
- のとき
  - ここで $\frac{1}{2}$ がでてくるわけでキョリは関係ない？
  - 偏微分係数の計算にナナメのマスを使っているだけ
さらにについてはに限らず以下のように与える
- $\tan\theta = \frac{\Delta y}{\Delta x}$
その上で
- $\phi = \frac{\pi}{2} - \theta$ として
- $\sin \phi \cos \phi$ をもう片方の係数とするべきだろう
ここで
- $\frac{\par^2}{\par y \par x}$ と $\frac{\par^2}{\par x \par y}$ が（ななめの）どちらの方向を表すのか分かっていないので、シミュレーションでたしかめられるかどうか
- $\tan\theta = t$ が与えられれば $\sin \theta \cos \theta =\frac{t}{1+t^2}$ 、 $\sin \phi \cos \phi =\frac{-t}{1+t^2}$ で計算できる