数値計算@常微分方程式@その2@オイラー法

参考文献及びWEB

参考文献

(1)「電気・電子工学のための数値計算法入門」橋本修著　総合電子出版社
(2)「ディジタル信号処理技術」玉井徳みち、長島厚、藤田泰弘、若井修造著日経BP社
(3)「よくわかる有限要素法」福森栄次著 Ohmsha

参考WEB

http://szksrv.isc.chubu.ac.jp/java/physics/rlc/rlc0.html

http://www.asahi-net.or.jp/~jk2m-mrt/kiso_RLC.htm

http://ja.wikipedia.org/wiki/%E4%B8%89%E8%A7%92%E9%96%A2%E6%95%B0

http://www.cmplx.cse.nagoya-u.ac.jp/~furuhashi/education/CircuitMaker/chap1.pdf

http://www.akita-nct.ac.jp/~yamamoto/lecture/2003/5E/lecture_5E/diff_eq/node2.html

http://chemeng.on.coocan.jp/cemath/cemath08.html

http://homepage3.nifty.com/skomo/f6/hp6_3.htm

http://homepage1.nifty.com/gfk/rungekutta.htm

あわせて読む

(a) 微分の説明

微分方程式の解を求める

(b) 微分方程式の解を数値計算で求める

積分@その1@積分とはなにか

(d) 関数が微分の形 f'(x) の形で表されているときの積分値を求める

(一般的に数値積分をする意味は、関数自体が変化の式(例えば速度の関数)などで表されているとき
ある時刻から時刻までの距離を求めたいときなどで使用される。

もちろんその変化の式が、手計算によって計算可能ならば積分によって時間と距離の関係を導出し
その時間を入れることで直接距離を導いたほうが早い。

しかし、一般的に時間と距離の関係のような式自体が導出することが不可能な場合が世の中には多く
時間と速度の関係のような変化の式を見つけ出すことの方が容易なのだが、これを手計算で積分すること
は非常に難しいため数値積分が提案された)

はじめに

微分方程式の初期解を求める方法として、オイラー法を説明する。

1階の常微分方程式は、一般に

$\frac{dy}{dx} = f(x,y)$ 　(1-1)

という形で表すことができる(微分が $\frac{dy}{dx}$ が、 f(x,y) と同じ意味になるのだが、
違和感があるが仕方がない!)。

いま、と x+h の区間内でテイラー展開すると、先に示したように $h \lt 1$ の範囲で

y(x + h) $= f(x) + h f'(x) + h^2 \frac{1}{2!}f''(x)+h^3 \frac{1}{3!}f'''(x) + \cdots$ (1-2)

ここで、を十分小さくすると

f(x+h) = f(x) + hf'(x) (1-3)

となり、

$y_{n+1} = y_{n} + h f(x_n,y_n)$ (1-4)

を得る。

ここで、初期条件 x = x_0 のとき y = f(x_0) = y_0 とすると、
(1-4)式から

$y_{1} = y_{0} + h f(x_0,y_0)$ (1-5)

次に x = x_1 = x_0 + h とし、それに対して y = f(x_1) = y_1 とすると

$y_{2} = y_{1} + h f(x_1,y_1)$ (1-6)

この考えを繰り返すと

$y_{3} = y_{2} + h f(x_2,y_2)$
$y_{4} = y_{3} + h f(x_3,y_3)$
$y_{5} = y_{4} + h f(x_4,y_4)$

$y_{n+1} = y_{n} + h f(x_n ,y_n)$ (1-7)

この刻み幅を極めて小さくすると解の精度は向上する。

このようにテイラー展開の1次の項のみ使用して解を求めるとき
1次の精度をもつアルゴリズムと呼ぶ。

修正オイラー法

増分を n, n+1 での微分の平均を取ることで表してみる

$y_{n+1} = y_{n} + \frac{f(x_n,y_n) +f(x_{n+1},y_{n+1})}{2} h$ (1-8)

たとえば、 x=x_0 , y=y_0 から始まる場合

$y_1 = y_0 + \frac{f(x_0,y_0) +f(x_{1},y_{1})}{2} h$

この時点で、右辺の y_1 は求められていないので、
この場合は、

y_1 = y_0 + h f(x_0,y_0)
$y_1 = y_0 + \frac{f(x_0,y_0) +f(x_{1},y_{0} + hf(x_0,y_0))}{2} h$

この考え方を一般形に変形すると、

$y_{n+1} = y_0 + \frac{f(x_0,y_0) +f(x_{1},y_{0} + hf(x_0,y_0))}{2} h$

さらにこれを簡単に表すと、

$y_{n+1} = y_n + \frac{y'_n + y'_{n+1}}{2} h$

オイラー法の例題

(1)SW1を入れたときの微分方程式を考える。

$L\frac{di}{dt} + Ri = E$ (1-1)

この微分方程式の一般解は、

特解(定常解@過度状態が終わった状態)のときの

i_s (1-2)

と
式(1)で E=0 とおいた場合の一般解(過渡解@過度状態のときの解)

i_t (1-3)

の和で表される。

i = i_s + i_t (1-4)

定常解(過度状態が終わった状態)を求める

定常解 i_s を以下のようにしてもとめる。

は定数であるため、 i_s は定数となる。

そのため時間変化がないから
$\frac{di_s}{dt}=0$
となるため(1)は、

Ri_s = E
$i_s = \frac{E}{R}$ (1-5)

過渡解(過渡状態)を求める

次に過渡解 i_t は、

$L\frac{di_t}{dt} + Ri_t = 0$
$L\frac{di_t}{dt} = - Ri_t$
$\frac{1}{i_t} di_t = -\frac{R}{L} dt$

両辺を積分すると

$\int \frac{1}{i_t} di_t = - \int \frac{R}{L}dt + \log A$
$\log i_t = -\frac{R}{L} t + \log A$

$i_t = Ae^{-\frac{R}{L} t}$ (1-6)

一般解

$i = i_s + i_t = \frac{E}{R} + Ae^{-\frac{R}{L} t}$

初期条件は、t=0のとき i=0 なので、

$0 = \frac{E}{R} + A$
ゆえに

$A = -\frac{E}{R}$

つまり

$i(t) = \frac{E}{R}(1-e^{-\frac{R}{L} t} )$ (1-7)

時定数 $\tau = \frac{L}{R}$ を代入して

$i(t) = \frac{E}{R}(1-e^{-t/\tau} )$ (1-7b)

この式は時間の経過とともに過度状態が終わり電流は、
$\frac{E}{R}$ に収束することを示している。

(2) SW2に切替えたときの方程式を求める

$L\frac{di}{dt} + Ri = 0$

$i = Ae^{-\frac{R}{L} t}$ (1-8)

(a) のとき

(1-7)式に t=s を代入したものと(1-8)式に t=s を代入したものは等しいので

$\frac{E}{R}(1-e^{-\frac{R}{L} s} ) = Ae^{-\frac{R}{L} s}$

$\frac{E}{R}(1-e^{-\frac{R}{L} s} ) * (e^{\frac{R}{L} s}) = Ae^{-\frac{R}{L} s} ) * (e^{\frac{R}{L} s})$

$A = \frac{E}{R}(e^{\frac{R}{L} s}-1 )$ (1-9)

となる。

ゆえに (1-8)式は、

$i = \frac{E}{R}(e^{\frac{R}{L} s}-1)*e^{-\frac{R}{L} t}$

$i = \frac{E}{R}(e^{\frac{R}{L} (s-t)}-e^{-\frac{R}{L} t} )$

$i = \frac{E}{R}(e^{-\frac{R}{L} (t-s)}-e^{-\frac{R}{L} t} )$

時定数 $\tau = \frac{L}{R}$ を代入して

$i = \frac{E}{R}(e^{-(t-s)/\tau}-e^{t/\tau} )$ (1-10)

となる。

プログラム

(1)SW1を入れたときのオイラー法

$L\frac{di}{dt} + Ri = E$
$L\frac{di}{dt} = E -Ri$
$\frac{di}{dt} = \frac{E-Ri}{L}$

$i_{n+1} = i_{n} + \Delta{h}\frac{E-Ri_{n}}{L}$

(2)SW2を入れたときのオイラー法

$L\frac{di}{dt} + Ri = 0$
$L\frac{di}{dt} = -Ri$
$\frac{di}{dt} = \frac{-Ri}{L}$

$i_{n+1} = i_{n} + \Delta{h}\frac{-Ri_{n}}{L}$

(注意)
(1)から(2)に以降する時間は、
(1)で十分に安定してからなので

if( t <= sw )

の状態を持つ。

コード

#include "stdafx.h"
#include<iostream>
#include<fstream>
#include<io.h>
#include<math.h>

double f(double x);
double kai(int in);

int _tmain(int argc, _TCHAR* argv[])
{

	errno_t err; 
	FILE *fp_output;
	
	double I=0;                    //電流
	const static double E=1.0;     //電圧源   
	const static double R=1.0;     //抵抗
	const static double L=1.0;     //コイル
	const static double tau = L/R; //時定数
	const static double dh=1.0e-1; //微小時間0.1
	const static double sw=1.0;    //両スイッチ切換え時刻

	//変数
	double t=0;    //実時間
	double i1=0.0; //初期電流(厳密解)
	double i2=0.0; //初期電流(オイラー法)
	double i3=0.0; //初期電流(修正オイラー法)


	if( err = fopen_s(&fp_output,"data.txt","w")  != 0){ exit(2);}
	
	
	for(int n=1; n< 5*(int)(sw/dh);n++){
		
		t=(double)n*dh;
		// 0 <= t <= s の場合
		if( t <= sw ){
			i1=(E/R)*(1.0-exp(-t/tau));
			i2=i2+dh*((E-R*i2)/L);
			i3=i3+dh*(((E-R*i3)/L)+((E-R*i2)/L))/2.0;
		}else{
			i1=(E/R)*( exp(-(t-sw)/tau) - exp(-t/tau) );
			i2=i2+dh*(-R*i2/L);
			i3=i3+dh*((-R*i3/L)+(-R*i2/L))/2.0;
		}

		printf("t = %f i1 = %f i2 = %f i3 = %f\n",t,i1,i2,i3);
		fprintf_s(fp_output,"%f %f %f %f\n",t,i1,i2,i3);

	}

	if(fp_output != NULL) fclose(fp_output);
	getchar();
	return 0;
}

結果の描画

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5.2 オイラー法
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'data.txt'

t=data(:,1);
i1=data(:,2);
i2=data(:,3);
i3=data(:,4);

plot(t,i1,'r-','linewidth',2);
hold on
plot(t,i2,'b+-','linewidth',2);
hold on
plot(t,i3,'g+-','linewidth',2);
hold on
grid on

xlabel('時間(sec)','Fontsize',18)
ylabel('電流(A)','Fontsize',18)

h=legend('厳密解',...,
'オイラー法',...,
'修正オイラー法',...,
1);
set(h,'FontSize',15)

axis([0 5.0 -Inf 0.8]);
h=gca
get(h)
set(h,'LineWidth',2,...,
'FontSize',18)

print -djpeg euler_method.jpg