微分方程式モデリングとシミュレーション

(1)

微分方程式

モデリングとシミュレーション

2018年度

1

(2)

質点の運動のモデル化



粒子と粒子に働く力



粒子の運動→粒子の位置の時間変化

粒子の位置の変化の割合→速度 速度の変化の割合→加速度



力と加速度の結び付け

Newtonの運動方程式：微分方程式 解は、時間の関数としての位置 ©只木進一（佐賀大学） 2

(3)

Newtonの運動方程式



質点の運動は、Newtonの運動方程式

で記述される

加速度は力に比例する ©只木進一（佐賀大学） 3 2 2 d d x m t = F  _

(4)

微分



独立変数𝑡𝑡とその従属変数𝑥𝑥(𝑡𝑡)



一階微分

変数𝑡𝑡が微小に増加した際の、𝑥𝑥の増分の割合：一定とは限らない 関数𝑥𝑥(𝑡𝑡)の時間に関する変化率 _{(𝑡𝑡, 𝑥𝑥)空間内の関数𝑥𝑥(𝑡𝑡)の傾き} 4 0 d ( ) ( ) lim d h x x t h x t t → h + − = ©只木進一（佐賀大学）

(5)



二階微分

関数𝑥𝑥(𝑡𝑡)の変化率𝑣𝑣(𝑡𝑡)の変化 関数𝑥𝑥(𝑡𝑡)の曲り具合（上に凸、下に凸など） 5 0 d ( ) d d ( ) ( ) lim d h x v t t v v t h v t t → h = + − = ©只木進一（佐賀大学）

(6)

微分方程式



関数の変化を記述することで、それに

従う関数を求める



例1：変化の割合が一定



例2：変化の割合が関数自身に依存

6 d d x a t = x t( ) = at +b d d x ax t = x t( ) = bexp

( )

at 微分方程式で定まらない定数（積分定数）𝑏𝑏が現れる ©只木進一（佐賀大学）

(7)

二階微分方程式



独立変数による二階微分を含む方程式



例：

7 2 2 2 d d x x t = −ω ( ) cos( ) sin( ) x t = a ωt +b ωt 微分方程式で定まらない定数（積分定数） 𝑎𝑎と𝑏𝑏が現れる ©只木進一（佐賀大学）

(8)

二階微分方程式を一階連立微分方

程式へ



_{𝑣𝑣 𝑡𝑡 = d𝑥𝑥(𝑡𝑡)/d𝑡𝑡を導入}

8 2 2 2 d d x x t = −ω 2 d d d d v x t x v t ω = − = ©只木進一（佐賀大学）位置𝑥𝑥の時間変化は 𝑣𝑣で定まる

(9)

数値積分の観点で微分方程式を

見ると



微分方程式を解くことを「積分」と呼

ぶ



微分方程式は、「左辺(導関数)が右辺

で与えられる」と読む

導関数→微小変化量→直後の関数の変化値 ©只木進一（佐賀大学） 9

( )

d , d y t t = f y  _ _

(10)

微分方程式を数値的に解くとい

うこと



微分方程式



ある𝑡𝑡での𝑥𝑥(𝑡𝑡)と𝑓𝑓(𝑥𝑥, 𝑡𝑡)が分かれば、

𝑥𝑥(𝑥𝑥 + Δ𝑡𝑡)の値を近似的𝑂𝑂( Δ𝑡𝑡

2

)に得る

ことができる

10

( )

d , d x f x t t =

(

)

( )

2 2 d ( ) d ( ) , x x t x t t t t O t t O t x t f x t ∆ ∆ + ∆ + ∆ = + + ∆ + = ©只木進一（佐賀大学）

(11)

Euler法



最も簡単な微分方程式の数値解法



一元一階微分方程式の場合



近似的な時間発展式

11

( )

d , d y f t y t =

(

t h

)

y t

( )

(

t,

( )

t

)

y +  + ×h f y ©只木進一（佐賀大学）

(12)

Euler法



_{𝑛𝑛元一階微分方程式の場合}



近似的な時間発展式

右辺は時刻𝑡𝑡における量だけで表現されている 12

( )

d , d i i y f t y t = 

(

)

( )

(

,

( )

)

i t h y ti h f t yi y +  + ×  t ©只木進一（佐賀大学）

(13)

Javaを用いた微分方程式の数

値解法



数値解法：Euler法またはRunge-Kutta法

微分方程式に依存しない一般的手法



微分方程式をインターフェイスのイン

スタンスとして渡す

C/C++の関数ポインタに相当 ©只木進一（佐賀大学） 13

(14)

微分方程式を表すインターフェ

イス



独立変数𝑡𝑡と従属変数𝑦𝑦

_𝑖𝑖

_𝑡𝑡

引数の値から、導関数の値を返す関数 14 @FunctionalInterface

Public interface DifferentialEquation {

public double[] rhs(double t, double y[]); }

©只木進一（佐賀大学）

関数インターフェイスであること示す注釈

(15)

Euler法によって数値積分を行

うクラス



微分方程式を引数で渡す

_{ℎだけ独立変数を進める}

15

public class Euler {

public static double[] euler(double t,double y[], double h, DifferentialEquation eq){

int n = y.length;

double yt[]=new double[n];

double dy[] = eq.rhs(t, y); for(int i=0; i<n ; i++){

yt[i] = y[i] + h * dy[i]; } return yt; } ©只木進一（佐賀大学） ( ) ( )

(

, ( )

)

i t h y ti h f t yi y + = + ×  t 微分方程式独立変数𝑡𝑡と従属変数 ⃗𝑦𝑦 からd ⃗𝑦𝑦/d𝑡𝑡を求める

(16)

Javaで微分方程式を定義する



DifferentialEquationはInterface

インスタンスを生成できない



二つの方法

(17)

関数の定義：一様重力の例

抽象クラスの実装を使って

public double[] rhs(double t, double[] y){ double dy = new double[n];

dy[0] = y[1];//dx/dt dy[1] = g; // dv/dt return dy; }; d d d d x v t v g t = =

(18)

関数の定義：一様重力の例

ラムダ式を使うと

= (double t, double y[]) -> { double dy = new double[n]; dy[0] = y[1];//dx/dt dy[1] = g; // dv/dt return dy; }; d d d d x v t v g t = = 各微分方程式をどのインデクスに割り当てるかは定めはない

(19)

Javaにおけるλ式



関数を表すインターフェースの実装を

簡潔に表す方法

関数そのもの（後で引数に値が入り評価される）を表現する 関数を変数として扱える

(引数並び) -> {関数の実体};



Java8以降で利用可能になった

(20)

20

import java.util.function. DoubleFunction; public class LambdaMain {

public static double generalSum(

double data[], DoubleFunction<Double> op){ double sum=0.; for(double d:data){ sum += op.apply(d); } return sum; }

public static void main(String[] args) { double data[]={1.,5.,8.,11.};

//二乗の和を計算

double sum = generalSum(data, d->d*d); System.out.println(sum);

} }

データに対する具体的演算は未定義

(21)

基本的な関数の定義例

java.util.functionに定義済み



DoubleFunction<R>

double型の一変数に対してR型の値を返す



UnaryOperator<T>

T型の一変数に対してT型の値を返す



Function<T,R>

T型の一変数に対してR型の値を返す



Predicate<T>

(22)

拡張されたfor ループ

(23)

ラムダ式の活用

Listなどの処理

©只木進一（佐賀大学） 23 List<Integer> list; // 要素を印刷 list.stream().forEachOrdered( x -> { int i = list.indexOf(x); System.out.println(i + ":" + x); } );

(24)

24

List<Integer> list; // 要素の和

int sum = list.stream().reduce(0,

(acc,item)->acc+item); // 要素の最大値

int max = list.stream().reduce(list.get(0),

(acc,item)->Math.max(acc,item); // 正の要素の数

微分方程式 モデリングとシミュレーション