2005/11/22 数値解析学と演習 福井 先生 今週と来週で常微分方程式をやる 今週 : オイラー法 次週 : ルンゲ=クッタ 常微分方程式 dy/dx = f(x,y) の形 # 前回の、数値積分と違いは、右側の式が、 f(y) から f(x,y) になったこと 解は、 y = \int f(x,y) dx + C の形 ( 積分定数 C が含まれる )。 常微分方程式には、一般に、積分定数が出てくるので、解が一意に定まらない => 数値的に解くには、積分定数の値を决めるための条件が必要 一階の常微分方程式 (積分定数が一つ) 1 点 ( 最初の点 ) を决める => 初期値問題 [例] 宇宙船の軌道を决める ニ階の常微分方程式 (積分定数が二つ) 2 点 ( 始点と終点 ) を决める => 境界値問題 [例] 電線の張り 両側 ( 境界 ) に重りを付ける # 線自身も重さを持つ # (電線の素材としての)銅とアルミ # 銅の方が伝導率は高いがアルミの方が軽い # 同じ太さなら銅が良い # 同じ重さならアルミが良い オイラー法/ルンゲ=クッタ => 初期値問題 境界値問題の解法は、偏微分方程式の解法に近い # 境界値問題を「初期値問題+α(例えば、二分法)」で解くことができる # 大砲を射って、はずれたら、補正する # より一般には.. 複数の条件を連立一次常微分方程式を立て解く オイラー法/ルンゲ=クッタ 一階の微分方程式を解くために設計されている 一階の微分方程式の解 発散するか減衰するかどちらか ( 一種の指数関数 ) 例 : 放射能の半減期 # y = e^x 二階の微分方程式の解 振動解になる # y = \sin{x}, \cos{x} # 実数解だと分かれるが、複数素数解としては同じ ## cf. e^{x+yi} = e^x ( \cos{y}+\sin{y}i ) == 一般の二階の微分方程式の解法 一般の形 d^2y/dx^2 + α dx/dt + βx = 0 # x は長さで、t は時間だと思っている 正規形に変換する ( cf. text p.472 ) x : 長さ t : 時間 dx/dt = v : 速度 d^2x/dt^2 = (d/dt)(dx/dt)=(d/dt)v: 可速度 この v を利用して、元の式を書き換えると、次のような連立微分方程式になる。 dx/dt = v dv/dt = -αv - βx # 形式的に、一階の微分方程式の形にできる # ただし、x の初期値の他に v の値の初期値も必要になることに注意 # 同様に三階微分でも三元連立の形にできる ( まあ、普通三階の微分方程式というものはないが.. ) この二階の微分方程式の解は、ほぼ、減衰 ( あるいは発散 ) する、振動解 近似解では、誤差が含まれる 振幅の誤差 : a 位相の誤差 : b # 位相のずれ 正弦波の場合 a\sin{t+b} a, b に誤差が含まれる 位相を合せようとすると ( 誤差があるので.. ) 振幅にしわ寄せが行く 振幅を合せようとすると ( 誤差があるので.. ) 位相にしわ寄せが行く 一般の近似解法では、両方に出るのが普通 位相の方が重要な場合 ( 建築関係の場合 ) => ニュマークのβ法 位相のずれをなくす工夫がされている # 結果的に振幅の誤差が大きくなる == [Overview] 常微分方程式 初期値問題 一階用の解法 オイラー/ルンゲクッタ # 二階以上は、正規形にすると、一階にできる 本質的には、指数関数の近似 # 複素数まで広げると、振動解も含まれている 近似法によって、振幅か位相 ( あるいは、その双方 ) に誤差 境界値問題 == [オイラー法] dy/dx = f(x,y) # 左側を差分化 ( 一階の近似式 ) \frac{y(x+h)-y(x)}{h) = f(x,y) # 整理 y(x+h) = y(x) + h f(x,y) # とくに x = 0 の時 y(h) = y(0) + h f(0,y(0)) # y(0) が初期値 y(0) がきまれば、これから y(h) が決る。 y(h) がきまれば、これから y(2h) が決る ... # これがオイラー法 # テーラー展開と比較すると論理的な誤差 ( 打ち切り誤差 ) は二次の以降項目の部分 y(x+h) = y(x) + \frac{h}{1!}y'(x) + \frac{h^2}{2!}y''(x) + + \frac{h^3}{3!}y^(3)(x) + .. # ほぼ、\frac{h^2}{2!}f''(x) の項目 ## O(h^2) の誤差 ### h を小くすれば、誤差が小くなる 刻み ( h ) を小くすると、その区間での誤差は確かに小くなる ( O(h^2) ) しかし、h を小くすると、繰り返し回数が多くなるので、誤差が大きくなる。 結局、区間全体では、誤差は O(h) になる => 一次式の近似なので、誤差が O(h) なるのは当然 # 区間 [a,b] での誤差評価 : e(a,b) \le c h ( c は関数によって異る定数値 ) ## 誤差評価の場合はどうしても、「良くわからない定数」(上記の c)が出てくる ### しかし、いずれにせよ h に無関係な値であることが重要 近似の次数を上げるべきか ? 解が、非常に滑らかであることがわかっていれば、次数を上げることに意味がある => 宇宙船の軌道計算 (ブッチャー博士) 八次, 九次の公式を作っている その代わりに、計算を速くするために刻みを大きくする # 宇宙の外に出てしまえば、純粋に、重力の影響のみ ## 打ち上げの場合は、空気の擾乱があるので、刻みを大きくとれないので、低次の公式を使わざるを得ない [数値解法の次元] オイラー法の次数 : 1 改良すると 2 次式まで ルンゲ=クッタの場合は 4 次式になる # 常微分方程式の誤差解析の問題点 # 誤差が二箇所に影響する # y(x+h) = y(x) + h f(x,y(x)) # ^^^^ ^^^^ ## 数値積分に比べて、誤差解析が大変 == 常微分方程式の良いところ 正規形にして、連立にすると微分の出る所と関数の所が分離できる 一般化が行いやすい == [今週の演習] 番号 11/22 問題 常微分方程式 y' = 2 x + 3 y y(0)=1 を オイラー法 で解く 区間 [0,1] 差分 h = 0.1 答(例) Δx = 0.1 の時 y(0.1) = 1.3 y(0.2) = 1.71 y(0.3) = 2.263 ... となる。 厳密解 y = \frac{11}{9}\exp{3x}-\frac{2}{3}x-\frac{1}{9} # y(0.3) の値を比較すると、12% 位の誤差ができる == 修正オイラー法 => 中心差分を使う 二次の近似になる 日本の常微分方程式の大家 田中正次 先生( ずっとルンゲ=クッタの色々な公式を作った ) # 以前、日大の工学部 == オイラー法の公式 y(x+h) = y(x) + h f(x,y) これは、銀行の利子の計算式 ( 複利計算 ) => 刻みを小くすると額が大きくなる # 真値を下から近似している