オイラー法による微分方程式の数値解法

dy/dx=-2xy^2
y(0)=1

でx=1での値の近似値をオイラーの方法で、求めよ(n=10)

という問題ですが、ウェブでオイラーの方法についてあらかた調べたのですが、記述が複雑うまく理解できませんでした。まず増やしていく幅のhは自由に設定していいのでしょうか?

オイラーの方法による解き方をやさしく教えていただけたら嬉しいです。

よろしくお願いします。

No.1 回答者： Tacosan 回答日時： 2011/06/09 23:22

問題の最後にある「(n=10)」の n って何?

No.2 回答者： muturajcp 回答日時： 2011/06/10 03:54

dy/dx=-2xy^2

∫y^{-2}dy=∫-2xdx
-y^{-1}=-x^2+c
-1=c
-y^{-1}=-x^2-1
y=1/(x^2+1)
y(1)=1/2

No.3 ベストアンサー 回答者： noname#154783 回答日時： 2011/06/10 04:58

1階の常微分方程式

dy/dx = f(x, y)
の解
y = y(x)
に対して，初期条件
y0 = y(x0)
が与えられているなら，十分小さいhに対して，
y(x0 + h) ≒ y(x0) + y'(x0) h = y0 + f(x0, y0) h
が成り立ちます（1次近似）．

# この1次近似がオイラー法の根本的な原理であり，これが解らなければ，オイラー法が解っていないというよりも，微分という考え方がまだ理解できていないと思われますので，もしそうなら，1次近似について復習してください．

そこで，
x1 = x0 + h
と置くと，
y(x1) ≒ y0 + f(x0, y0) h = y1
と表せます．

この式で求められるy1は飽くまで近似値なので，真の解のグラフが点(x1, y1)を通るとは限らないのですが，当たらずとも遠からずってとこでしょうから，点(x1, y1)が真の解のグラフ上にあるものとみなし，この点で同じ近似を行います：
y(x2) ≒ y(x1) + y'(x1) h = y1 + f(x1, y1) h = y2

同じようにして，逐次近似を行っていくと，一連の近似値が得られます：
y1 = y0 + f(x0, y0) h,
y2 = y1 + f(x1, y1) h,
y3 = y2 + f(x2, y2) h,
y4 = y3 + f(x3, y3) h,
...

そこで，xにおける解の値y(x)の近似値が欲しければ，区間[x0, x]をn等分し（本当は等分でなくてもいいのですが，簡単のためそうします），それぞれの分点に
x0, x1, x2, ..., x[n-1], xn = x
と名前を付けると，
y1 = y0 + f(x0, y0) h,
y2 = y1 + f(x1, y1) h,
...
y[n-1] = y[n-2] + f(x[n-2], y[n-2]) h,
yn = y[n-1] + f(x[n-1], y[n-1]) h.

これがオイラー法です．

このやり方は素朴で分かりやすいのですが，誤差が蓄積しますので，精度はあまり良くありません．

で，今回の微分方程式
dy/dx = -2x y^2,
y(0) = 1
ですが，VBSで簡単なコードを書いてみました：


'Euler法

Option explicit

Function f(x, y)
f = -2*x*y^2
End Function

Dim x, xn, y, h
Dim n, i

'初期条件
x = 0
y = 1

xn = 1 '評価点

n = 10 '分割数
h = (xn - x)/n '幅

For i = 1 To n
y = y + h*f(x, y)
x = x + h
Next

MsgBox "y(" & xn & ") = " & y


これをメモ帳でも何でもいいからテキストエディタで拡張子.vbsのテキトーなファイル名で保存し，アイコンをダブルクリックすると，

y(1) = 0.503641976039014

とか値が出力されます．

この微分方程式は，ANo.2さんが回答してくださっているように解析的に解けるのですが（変数分離形），真の解は
y(1) = 0.5
なので，確かに「当たらずとも遠からず」って感じです．