オイラー法は微分方程式を解く

3538 ワード

by Conmajia 2014
原文は私が2014年に書いたもので、C#で完成しました.ここではJavaScriptに変更しました.
特基礎.もちろん最も便利なのは、数学ライブラリやMatlab、Mathematicsなどの数学ソフトウェア(値だけを求めている場合)、C、Java、Go、Erlangの他の言語に変えて実現することです.Euler(Euler)と中心差分が近づくのは,最も素朴な考え方であり,残念ながら代数精度が低すぎるが,竜格庫塔の安定性はまた問題である.要するに普通のものを計算するのにしか使えない.高大上の問題は一般的に一般化$alpha$,Wilson-$Theta$などを用いて、係数を利用して人工減衰を調節し、低周波の部分をアルゴリズムの影響を少なくし、高周波減衰を増大させ、微分方程式の安定性を増加させ、音響アルゴリズムの精度を低下させる.
Math.NET計画は.NET用の数学演算ライブラリであり、数値計算ライブラリIridium.NETを含む.これは
非常に古いライブラリで、最近の更新は10年前です.でも、
数学のものは永遠に時代遅れではないので、よく研究してみると、干物がたくさん食べられます.$ git clone git://github.com/mathnet/mathnet-iridium.gitこのような微分方程式を解く必要がある場合があります.
\[ y'=y-\cfrac{2x}{y}\]
(y(0)=1).両側積分,得られる:
\[\begin{align*}\int_{x_n}^{x_{n+1}}y'\mathrm{d} x&=\int_{x_n}^{x_{n+1}}\left(y-\cfrac{2x}{y}\right)\mathrm{d} x\\y_{n+1}-y_n&=\int_{x_n}^{x_{n+1}}\left(y-\cfrac{2x}{y}\right)\mathrm{d}x\\y&=\sqrt{1+2x}\end{align*}\]
工程の上でただ数値解だけが要って、一般的に微分の代わりに差分近似を採用します.最も簡単で最も素朴な方法はオラの公式を使うことです:
\[\tag{1} y_{n+1}=y_n+\Delta x f(x_n,y_n)\quad n=0,1,2,\cdots\]
導出は簡単です(x_n)、あります:
\[ y'_n=f(x_n,y_n).\]
はい(y'(x_n))あります.
\[\tag{2} y'_n\approx\cfrac{y_{n+1}-y_n}{\Delta x}.\]
式(2)左は微商(喜提迪丽热巴のあの微商ではありません)、右は差商(淘宝伪劣制品の売り手ではありません)、(Delta x=left|x_{n+1}-x_nright|).
代回式(2)、あります:
\[\begin{align*} y'_n=f(x_n,y_n)\Rightarrow\cfrac{y_{n+1}-y_n}{\Delta x} &\approx f(x_n,y_n)\\y_{n+1}-y_n &\approx\Delta x f(x_n,y_n)\\y_{n+1} &\approx y(x_n)+\Delta x f(x_n,y_n).\end{align*}\]
式(1)のオーラ式は次のようになります.
\[\tag{3} y_{n+1}=y_n+\Delta x\left(y_n-\cfrac{2x_n}{y_n}\right).\]
残りをmyFn関数で表すと仮定(y_n-dfrac{2 x_n}{y_n}):
function myFn(x, y) {
  return y - 2 * x / y;
}
Eulerは次のように実現できます.
function calculate(x0, y0, delta, xn) {
  var yn;
  while(x0 < xn) {
    yn = y0 + delta * myFn(x0, y0);
    y0 = yn;
    x0 = x0 + delta;
  }
  return yn;
}
今から実験を始めることができます.
理論上:
\[ y=\sqrt{1+2x}\Rightarrow y(1)=\sqrt{3}\approx 1.7321\]
(sqrt{3})は方程式の真の値であり、プログラムの目標は計算によって、できるだけ真の値に近い結果を得ることである.
(x_0=0),(y_0=1),(Delta x=0.1)をとり、プログラム計算結果は次のとおりです.
ans = 1.784771
誤差3.04%減少(Delta x)、たとえば(Delta x=0.0001)、計算結果は次のとおりです.
ans = 1.732112
誤差は0.0007%で、真値に非常に近い.
この方法は比較的正確な解を求めることができるが、(Delta x)が小さすぎるとwhileが何度も実行する、効率が低下する.
定積分の台形式を導入します.
\[\int_{x_n}^{x_{n+1}}f(x,y)\mathrm{d}x\approx\cfrac{\Delta x}{2}\left[f(x_n,y_n)+f(x_{n+1},y_{n+1})\right]\]
式(3)は次のように書くことができます.
\[\tag{4} y_{n+1}\approx y_n+\cfrac{\Delta x}{2}\left[f(x_n,y_n)+f(x_{n+1},y_{n+1})\right].\]
式(3)と式(4)の特徴は、前者は速度が速く、精度が低いことである.後者は速度が遅く、精度が高い.
\[ y'_{n+1}=y_n+\Delta x f(x_n,y_n)\]
精算、使用可能:
\[ y_{n+1}=y_n+\cfrac{\Delta x}{2}\left[f(x_n,y_n)+f(x_{n+1},y'_{n+1})\right].\]
これに合わせて、まず粗算により(y'{n+1})の近似値を得てから精算を行い、高精度の最終解を得る.これにより、計算結果の精度が保証されるとともに、計算資源があまり消費されず、計算効率が保証される.以下、改良された計算プログラムである.
function calculate(x0, y0, delta, xn) {
  var yp, yc;
  while(x0 < xn) {
    yp = y0 + delta * myFn(x0, y0);
    yc = y0 + delta * myFn(x0 + delta, yp);
    y0 = 1 / 2 * (yp + yc);
    x0 = x0 + delta;
  }
  return y0;
}
(x_0=0),(y_0=1),(Delta x=0.1)をとり、次のように計算されます.
ans = 1.737867
誤差0.33%
最も古いバージョン(誤差3.04%)と比較する、(Delta x)が同じである、サイクル数が近いことを意味する場合には、精度が10倍近く向上する.
具体的な问题の具体的な分析、手で微分方程式を求めるのは基本的にこのような考え方です.もちろん、「数値分析」と「工程数学」のこれらの课程にとって、私の上に书いたものは小児科にすぎません.もっと多くの时、手を伸ばす党を作るのは実はとてもいいです.
The End.\(\Box\)