• 締切済み

数値積分の誤差の問題

数値解析は初心者です。 惑星の運動を調べるために(3つの星です)、 微分方程式を数値的に解こうとしています。 エネルギーが保存する系では、 シンプレクティック法というもの方が良いと 教えてもらいいろいろな人のものを参考にして 使っています。 問題点は、シンプレクティック法の結果は、 ルンゲクッタよりエネルギーの保存の 精度はいいのですが 1~10万ステップくらいの挙動が 計算するコンピューターで違うことが 分かりました。誤差が出ているのだと思います。 ルンゲクッタは保存は悪い代わりに 10万ステップくらいまでマシンに依存した 結果はでません。 エネルギーは高い精度で保存しているので バグではないと思います。 シンプレクティック法は6次と8次を使っています。 どなたかシンプレクティック法を使って 不安定な系を扱っている方で シンプレクティック法の精度を上げる方法を ご存知の方がいらっしゃったらお教えください。

みんなの回答

  • tomtom_
  • ベストアンサー率39% (43/110)
回答No.1

情報が少ないので推測で書きます. まず,計算機によって結果が異なるというのが,非常に小さな変化なら気にすることはないと思います.シンプクレクティック数値解法や多様体法,シンメトリック法はそんなもんです. 次に,6次や8次という高い精度をどうやって作っているかが気になります.もしガウス・ルジャンドル型なら,反復回数を打ち切るマシンイプシロンなんかがどうなってるかチェックしたいところです. もし,対称分解法や対称合成法,フラクタル分解法を使っているなら,分解の仕方によって安定度が大きく変わります.例えば2次の解法から4次の解法を作る場合, S4=S2(a τ)S2(a τ) と分解するよりも, S4=S2(a τ)S2(b τ)S2(a τ) とか S4=S2(a τ)S2(b τ)S2(c τ)S2(a τ)S2(b τ) と分解するほうが,aやbなどの定数の値が小さくなって安定します.もちろん計算時間は増えます. もしこの辺りのせいなら,刻み幅の値を変えて計算し,比べてみれば分かります. 刻み幅を半分にして,誤差が6次の方法で0.5^6に減っていなかったら,プログラムのバグです. 誤差がちゃんと減っているば,丸め誤差や打切り誤差のせいですから,安定した解法を使うことは意味があります.

monohon
質問者

お礼

初心者にも分かりやすく、 詳しいご説明大変ありがとうございます。 誤差について詳しくないのでだいぶわかりました。 ありがとうございます。 そんなものと考えればいいとのことで 少し安心しました。 ただルンゲクッタより計算機の依存性が大きかったのが 意外でしたので。 私がもらったシンプレクティック法のやり方は、 6次と8次はすでに2次で分解した状態に解いた 時の各2時の係数の表が入っていてそれを 読み込んで使っています。これは、吉田という人が 論文で発表してあるものだそうです。 やはり時間ステップを小さくすれば誤差は減るんですか。今ステップを0.01単位時間にして、1万単位時間くらい計算してるのですが、ステップを0.001にした場合、1万時間まで見るためには、10倍のステップがいるために結局同じなのかなと勘違いしていました。 幅を変えるのをやってみます。 また報告します。 ありがとうございます。

monohon
質問者

補足

つかっているアルゴリズムは以下の通りです。 もっと良い分割の仕方や係数の表をご存じでしたらお教え頂ければと思います。 void symplectic_integrator_2(double *x, double dt) { double dxdt[2*DOF]; int i; /* Step 1 */ time_derivatives(x, dxdt); for (i = DOF; i < 2*DOF; i++) x[i] += 0.5*dxdt[i]*dt; /* Step 2 */ time_derivatives(x, dxdt); for (i = 0; i < DOF; i++) x[i] += dxdt[i]*dt; /* Step 3 */ time_derivatives(x, dxdt); for (i = DOF; i < 2*DOF; i++) x[i] += 0.5*dxdt[i]*dt; } void symplectic_integrator_6(double *x, double dt) { int i; double w0, w1, w2, w3; w1 = -1.17767998417887; w2 = 0.235573213359357; w3 = 0.784513610477560; w0 = 1.0 - 2.0 * (w1 + w2 + w3); /* Step 1 */ symplectic_integrator_2(x, w3 * dt); /* Step 2 */ symplectic_integrator_2(x, w2 * dt); /* Step 3 */ symplectic_integrator_2(x, w1 * dt); /* Step 4 */ symplectic_integrator_2(x, w0 * dt); /* Step 5 */ symplectic_integrator_2(x, w1 * dt); /* Step 6 */ symplectic_integrator_2(x, w2 * dt); /* Step 7 */ symplectic_integrator_2(x, w3 * dt); } void symplectic_integrator_8(double *x, double dt) { int i; double w0, w1, w2, w3, w4, w5, w6, w7; w1 = 0.311790812418427; w2 = -0.155946803821447*10.0; w3 = -0.167896928259640*10.0; w4 = 0.166335809963315*10.0; w5 = -0.106458714789183*10.0; w6 = 0.136934946416871*10.0; w7 = 0.629030650210433; w0 = 1. - 2. * (w1 + w2 + w3 + w4 + w5 + w6 + w7); symplectic_integrator_2(x, w7 * dt); symplectic_integrator_2(x, w6 * dt); symplectic_integrator_2(x, w5 * dt); symplectic_integrator_2(x, w4 * dt); symplectic_integrator_2(x, w3 * dt); symplectic_integrator_2(x, w2 * dt); symplectic_integrator_2(x, w1 * dt); symplectic_integrator_2(x, w0 * dt); symplectic_integrator_2(x, w1 * dt); symplectic_integrator_2(x, w2 * dt); symplectic_integrator_2(x, w3 * dt); symplectic_integrator_2(x, w4 * dt); symplectic_integrator_2(x, w5 * dt); symplectic_integrator_2(x, w6 * dt); symplectic_integrator_2(x, w7 * dt); }

関連するQ&A

  • 数値微分法についてです。

    数値的に微分を評価する時に、中心差分を使っているのですが、 どう考えても数値誤差としか思えない結果しかでません。 区切り幅は誤差が最も小さくなるように選んでいます。 中心差分よりも精度の良い数値微分法があれば教えていただけないでしょうか。 評価する関数は解析的に与えられておらず、補間して得られるようなものです。 (補間の精度にもよるのだと思いますが・・・)

  • オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法のC言語プログラムに

    オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法のC言語プログラムについて教えてください! 課題なのですが、まったくわからず困ってます>< 1 常微分方程式 dy/dx=f(x,y),y(0)=1 の数値解をオイラー法を用いて計算するプログラムを作為せよ。ただし、f(x,y)=3-6x^2-4x+2xyとする。 2 α=1,β=1,γ=1/2,σ=1/2 の場合の2次ルンゲクッタ法を考える。1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を求めるプログラムを作成せよ。また、オイラー法と2次ルンゲクッタ法の実行結果を示して、2つの近似精度を比較せよ。 3 1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を4次ルンゲクッタ法を使って求めるプログラムを作成せよ。また、オイラー法、2次ルンゲクッタ法、4次ルンゲクッタ法の実行結果を示して、3つの近似精度を比較せよ。 以上の3つです。 休んでいた自分が悪いのですが、ネットで調べてもよくわからなくて… わかる方、よろしくおねがいします…

  • 高階連立常微分方程式の数値計算

    4次のルンゲクッタ法を用いた数値計算を勉強しています. 1階連立常微分方程式と高階常微分方程式は理解でき,プログラムも作成することができました. 次に高階の連立常微分方程式を解こうと思ったら,頭が混乱してしまいました. 4次のルンゲクッタ法を用いて高解連立常微分方程式を解く考え方を教えて頂ければ嬉しいです. また何か良い参考書があれば教えて頂きたいと思います. よろしくお願いします.

  • 数値解析 微分方程式

    数値解析における常微分方程式を解くために用いる手法についてです。 オイラー公式、ホルン公式、ルンゲクッタ、ルンゲクッタ4次、有限差分法の関係 違いがよくわからなのでどなたか教えてください また常微分方程式を有限差分で解くとなったとき、結局オイラー公式などを使うと言う認識で間違いないでしょうか?それとも有限差分だけで解けるのでしょうか

  • 微分を数値的に積分して近似解を求める方法

    du/dt=f(u,t) というような感じの単純な常微分方程式を逐次的にコンピュータで数値積分していく場合、単純陽解法とかルンゲクッタ法とか予測子・修正子法、台形公式などいろいろあると思います。あるいは反復計算を含むものなどもあると思います。これらについてグラフを作図して解析解との差を評価したりするわけですが、偏微分方程式で使う場合も全く同じ評価がそのまま当てはまるでしょうか。常微分方程式での知見がそのまま偏微分方程式の計算に当てはまるのか、ということなのですが。反復計算を用いる方法で常微分方程式では芳しくない結果だったのですが、偏微分の場合は反復計算の価値が高くなるのではないかと思ったりするものですから。反復計算、すなわち繰り返しをやるうちに陰解法的になっていくからなのですが。 いかがでしょうか。細かい内容の質問で恐縮ですが、よろしくお願いします。

  • 減衰を持つ系のシンプレティック積分

     ず~っと、必ず減衰のある構造系分野に住んでいました。ある時、保存系の数値積分に迫られ、それに今まで日常的使っていた、ニューマークβ法やウィルソンのθ法を適用してみると、余りにあっけなく解が発散するので、びっくりしました。もちろんエネルギーを保存するように、β,γ,θを選びました。ルンゲ・クッタは多次元への拡張が面倒そうだったので、最初から無視していましたが、いずれそうなった(発散した)と思います。仕方ないので、小さな減衰を入れて対処しましたが、非常に不満が残りました。  その後、シンプレティック積分法というのを見つけ、やってみると、保存系に対して、ものすごい威力を発揮してくれました。以下簡単のため、一次元運動を考え、一般化運動量と一般化変位を(p,q)で表します。  シンプレティック積分法の解が、何故あれほど安定なのかを考えました。保存系であれば、エネルギー曲線E=H(p,q)があります。シンプレティック法は、一種の蛙跳び法ですが、pとqの各増分ステップにおいて、必ずエネルギー曲線Eを横切る事が、定性的に保証されているので、解が安定なのだと思います。そうすると同じ時間ステップ幅に対しては、ニューマークβ法やウィルソンθ法,ルンゲ・クッタ法よりも精度が落ちる事にはなりますが、発散しないという意味においては、シンプレティックの方が、よほど精度が良い事になります。  さらにシンプレティックは、アルゴリズムを簡単に出来ます。実際、保存系の一次のシンプレティック法は陽解法で、一次のオイラー法とアルゴリズムの手間は同等であり、簡単でありながら発散しないという、夢の陽解法です。精度を求めたいなら、時間Δtステップを狭くすればいいだけの話です。そして、シンプレティック法の適用は、読んだ限りでは保存系に限られていました。  しかし以上の状況を考えると、非保存系にシンプレティック法を直接適用しても、けっこういけるのでは?、と思えるのです。何故ならシンプレティック法の要は、各増分ステップにおいて、ΔpとΔqが、必ずエネルギー曲線を横切る事です。  十分短いΔtを取れば、「その瞬間のE」を「ΔpとΔq」が横切るように出来る、と思えます。駄目だったらΔtを短くします。この場合、計算ステップは増えますが、  ・シンプレティックは、アルゴリズムをオイラー法なみに簡単にできる.  ・現在のPCは、恐ろしく速くて、昔では考えられないくらいのメモリとディスクを搭載している. ことを考えると、自分には、十分魅力的な方法に映ります。この辺りの情報って、ないでしょうか?。ちなみに参考書(あまりないですが)は、  計算物理入門,臨時別冊・数理科学SGCライブラリ10,上田顕,サイエンス社,2001年. を使っています。  よろしくお願いします。

  • ルンゲクッタ法について

    ルンゲクッタ法についての問題なのですがよくわかりません。 常微分方程式 du/dt=u について、ルンゲクッタ法(3次.4次)による誤差分析を行えという問題です。

  • 精度と確率誤差

    10ml容ホールピペットの器差を求める実験を行なったんですが、この実験の精度を測ることができません。 各計算結果は以下のように出ました。 ホールピペットの20℃における内容積を調べたところ、9.974、9.973、9.980という結果が出た。 この平均値は9.976、標準偏差は14.11、確率誤差は9.517であった。 この中で上の確率誤差に着目して精度を測りたいんですが、イマイチ確率誤差の定義が分かりません。 確率誤差と精度はどのように関連していますか? この実験における確率誤差からは精度に関してどのようなことが言えますか? 有識者の方ご指南お願いします。

  • この方程式、stiffなんでしょうか・・・

    僕は以下のような連立微分方程式を数値計算で解こうと試みています。 現在微分方程式の解を級数展開と変数変換とルンゲクッタ法で求めようとしています。 しかし、どのように初期値をいじっても必ず解が発散してしまいます。 色々と自分で調べてみる過程で微分方程式を数値計算で解くときにstiff問題というものが生じることがあると知りました。 解の挙動が急に変わる場合や境界条件に大きく依存する場合などで途中からルンゲクッタ法などで数値計算できなくなり、解が発散してしまうという条件がぴったり自分の場合と一致しました。 以下の連立微分方程式はstiffなのでしょうか? (x-1)^4*f ''(x)+(2x-1)(x-1)^3/x*f '(x)-{ ((x-1)/x+A(x))^2+f^2-1 }f=0, (x-1)^4*A''(x)+(2x-1)(x-1)^3/x*A'(x)-(x-1)^2/x^2A(x)-f^2(A+(x-1)/x)=0 stiffだとルンゲクッタ法は使えなくなってしまうんですよね・・・・?

  • 測量での「誤差」と「精度」の検証について

    大学の「測量学」課題について教えてください 課題内容は 1.まず、ある五角形を、2つの方法で三角形に分割する   (1)「放射法」(三角形が5つできる)→【五角形1】とする   (2)五角形の1つの頂点から放射状に線を延ばして分割     (三角形が3つできる)→【五角形2】とする 2. 1で分割した五角形の面積を計算する   (1)【五角形1】を「ヘロンの公式」を用いて面積計算   (2)【五角形2】を「ヘロンの公式」を用いて面積計算   (3)【五角形2】を「S=(b×c×sinA)/2」(2つの辺と挟まれた角)     を用いて面積計算   結局、3つの面積の計算結果が得られる 3.上記1と2について、図上での測量及び面積の計算結果を比べ、   精度を検証して、その誤差について考察せよ 1、2は出来たのですが、3の意味がよくわかりません。 何か計算式を用いて考察するということでしょうか 図面はCADで作るので、元になる五角形は3つとも全く同じ大きさのものを使用しています 教科書には「○○誤差」という項目がたくさんあるのですが、どんな誤差が今回あてはまるのか、 「精度の検証」とはどうすればいいのか、 三角形の分割した時の辺と角についても誤差を考えよと言っているのか・・・。 実際に測量した数値の、誤差の処理法などはたくさん載っているのですが 異なる方法で計算した面積の考察は見当たらず、すごく困っています。 どのように考え、まとめるものなのか検討もつかない状態です ちなみに、測量学では小数点以下の処理に四捨五入は使わないものなのでしょうか? 通信教育なので授業がなく、聞く人がいない上に、レポート期限が迫っていてとても困っています 「恐らく○○を考えればいいんじゃ・・」みたいなことでもかまいませんので どなたかわかるかた、教えてください!!