5

私は Haskell を初めて使用し、演習として、Joel Franklin の著書Computational Methods in Physicsから (Mathematica で記述された) コードの一部を実装しようとしました。最初の引数としてラムダ式 (加速度) を受け取る次のコードを作成しました。一般に、加速度は x'' = f(x',x,t) の形式ですが、常に 3 つの変数すべてであるとは限りません。

-- Implementation 2.1: The Verlet Method
verlet _  _  _  _  0 = [] 
verlet ac x0 v0 dt n = take n $ zip [0,dt..] $ verlet' ac x0 v0 0 dt
  where verlet' ac x0 v0 t0 dt = 
          let
            xt = x0 + dt*v0 + 0.5*(dt^2)*(ac x0 v0 t0)
            vt = v0 + dt*(ac x0 v0 t0)
          in
            xt:(verlet' ac xt vt (t0+dt) dt)

ghci では、次のコマンドでこのコードを実行します (加速関数 a = -(2pi) 2 x は本の練習問題から来ています)。

verlet (\x v t -> -(2*pi)^2*x) 1 0 0.01 1000

私の問題は、これが実際には Verlet 法ではないことです。ここでは x n+1 = x n + dt*v n +1/2*a(x n ,v n ,n) ですが、Verlet 法はウィキペディアで説明されています。 x n+1 = 2*x n - x n-1 +a(x n ,v n ,n) になります。Verlet 積分法をより忠実に表現するには、この関数をどのように書き直せばよいでしょうか?

接線的に、これをよりエレガントかつ簡潔に書く方法はありますか? これを容易にする線形代数ライブラリはありますか? アドバイスありがとうございます。

4

3 に答える 3

5

忠実なベルレ列には、 x -- x n-1と x n-2の前の 2 つの値に応じて x nがあります。このような数列の標準的な例は、次のワンライナー Haskell 定義を持つフィボナッチ数列です。

fibs :: [Int]
fibs = 0 : 1 : zipWith (+) fibs     (tail fibs)
                        -- f_(n-1)  f_n

これは、フィボナッチ数列を無限 (遅延) リストとして定義します。への自己参照tail fibsは前の用語を示し、への参照はその前fibsの用語を示します。次に、用語を と組み合わせて(+)、シーケンス内の次の用語を生成します。

次のように、問題に対して同じアプローチを取ることができます。

type State = (Double, Double, Double)  -- (x, v, t) -- whatever you need here

step :: State -> State -> State
step s0 s1 = -- determine the next state based on the previous two states

verlet :: State -> State -> [State]
verlet s0 s1 = ss
  where ss = s0 : s1 : zipWith step ss (tail ss)

データ構造Stateは、必要な状態変数 (x、v、t、n、...) を保持します。この関数は、フィボナッチの場合にstep似ており、前の 2 つの状態から次の状態を計算します。(+)このverlet関数は、最初の 2 つの状態が与えられると、状態のシーケンス全体を決定します。

于 2014-03-09T06:50:12.770 に答える
0

user5402 からの提案を実装した後の私の解決策は次のとおりです。

-- 1-Dimensional Verlet Method
type State = (,,) Double Double Double -- x, x', t

first :: State -> Double
first (x, _, _) = x

second :: State -> Double
second (_, x, _) = x

third :: State -> Double
third (_, _, x) = x

verlet1 :: (State -> Double) -> State -> Double -> Int -> [State]
verlet1 f s0 dt n = take n ss
  where
    ss = s0 : s1 : zipWith step ss (tail ss)
      where 
        s1 = (first s0 + dt*(second s0) + 0.5*dt^2*(f s0), 
              second s0 + dt*(f s0), third s0 + dt)
        step :: State -> State -> State
        step s0 s1 = (2*(first s1) - first s0 + dt^2*(f s1), 
                      second s1 + dt*(f s1), third s1 + dt)

次のコマンドを使用して、ghci で実行しました。

verlet1 (\x -> -(2*pi)^2*(first x)) (1, 0, 0) 0.01 100

それはまさに私が期待していたものをもたらすようです-明らかに正弦波の動きです! 私はまだ x をプロットしていません (Haskell でそれを行う方法について何か提案があれば、歓迎します)。また、明らかなリファクタリングを見つけた場合は、遠慮なく指摘してください。ありがとう！

于 2014-03-10T07:04:43.130 に答える