数値計算@常微分方程式@その3@ルンゲクッタ法

参考文献及び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) の形で表されているときの積分値を求める

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

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

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

4次のルンゲクッタ法

(☆)注意微分方程式の解を求めるのであって、微分するのではない。
つまり、微係数を求めるのではない

ある連続関数 y=f(x) , $a\lt x \lt b$ において、区間[a,b]をn等分し、その間隔を $h=\frac{b-a}{n}$ とする。

$(\frac{dy}{dx})_i = \frac{ y_{i+1} - y_{i} }{h}$
の
$y_{i+1} = (\frac{dy}{dx})_i \times h + y_{i}$
の
$y_{i+1}$ を求めるのが常微分方程式の解を求める作業。

微分方程式の初期解を求める方法として、ルンゲクッタ法を説明する。

ルンゲクッタ法を簡単にいえば、適切な傾き(微係数)を求めて
(x_n , y_n ) と傾きおよび $x_{n+1}$ までの距離から $y_{n+1}$ を求めることである。

解法を以下に説明する。
(1) オレンジ色の直線の微分は、点 x_n , y_n の微分なので、

x_n での微係数は、

k1(n)=f(x_n , y_n ) (1-1)となる。

(2-a) 幅h/2の間、yの関数の傾きは, k1(n) であったと仮定して、 y(x_n+h/2) の値を仮定する。

$y(x_{n}+\frac{h}{2})=k1(n)*( \frac{h}{2}) + y_n$

そうすると、 $(x_n + \frac{h}{2} , y(x_n+\frac{h}{2}))$ の微係数は、

$k2(n)=f(x_n+\frac{h}{2},y_n+\frac{h}{2}*k1(n))$ (1-2)

(2-b) 次に再び幅h/2の間、yの関数の傾きは, k2(n) であったと仮定すると

$y(x_{n}+\frac{h}{2})=k2(n)*( \frac{h}{2}) + y_n$

そうすると、 $(x_n + \frac{h}{2} , y(x_n+\frac{h}{2}))$ の微係数は、

$k3(n)=f(x_n+\frac{h}{2},y_n+\frac{h}{2}*k2(n))$ (1-3)

次にこの傾きが十分に正しいと仮定して、実際に $x_{n+1}$ のときの $y_{n+1}$ を求めてみる。

$y(x_{n+1})=k3(n)*(h) + y_n$

そうすると、 $(x_{n+1} , y(x_n+h))$ の微係数は、

k4(n)=f(x_n+h,y_n + k3(n)*h) (1-4)

この4つの微係数に1:2:2:1の重みをかけて6で割った微係数にhをかけることで $y_{n+1}$ は求められる。

$y_{n+1} = y_n + \frac{1}{6}( k1(n) + 2k2(n) + 2k3(n) + k4(n) )*h$ (1-5)

まとめると以下となる。

k1(n)=f(x_n , y_n ) (1-1)
$k2(n)=f(x_n+\frac{h}{2},y_n+\frac{h}{2}*k1(n))$ (1-2)
$k3(n)=f(x_n+\frac{h}{2},y_n+\frac{h}{2}*k2(n))$ (1-3)
k4(n)=f(x_n+h,y_n + k3(n)*h) (1-4)

$y_{n+1} = y_n + \frac{1}{6}( k1(n) + 2k2(n) + 2k3(n) + k4(n) )*h$ (1-5)

4次のルンゲクッタ法のコード

■ $L\frac{di(t)}{dt} + Ri(t) = E\sin(\omega t)$ ときの解

は、

$\frac{di}{dt} + \frac{R}{L}i = \frac{E}{L}\sin(\omega t)$ (A-1)

となるので
$P(x) = \frac{R}{L}$
$Q(x) = \frac{E}{L}\sin(\omega t)$
n=0

となる。 $u = i^{(1-n)}$ とおいて変数変換すると

u' + (1-n)P(x)u = (1-n)Q(x)
$u' + \frac{R}{L}u = \frac{E}{L}\sin(\omega t)$ (A-2)

$\frac{E}{L}\sin(\omega t)$ が0のときの解をまず考える(過度解)

$\frac{du}{dt} + \frac{R}{L}u = 0$

$\frac{du}{u} = - \frac{R}{L}dt$
両辺を積分すれば、

$\int \frac{du}{u} = -\frac{R}{L}\int dt$

$\Longleftrightarrow \ln |u| = -\frac{R}{L}t + c$

$u = \pm e^{-\frac{R}{L}t+c}$

ここで $C = \pm e^{c}$ とすれば、
この方程式の解は

$u =C e^{-\frac{R}{L}t}$ (A-3)

( は0を含む任意の実数）となる。

■ 定数変化法で、(A-1)の一般解を得る

(A-3)の積分定数Cをxの関数とすると

$du = e^{-\frac{R}{L}t}{\rm dC} -\frac{R}{L}Ce^{-\frac{R}{L}t}{\rm dt}$ (A-4)

(A-4)に(A-2)(A-3)を代入する
$\frac{e^{-\frac{R}{L}t}{\rm dC} -\frac{R}{L}Ce^{-\frac{R}{L}t}{\rm dt}}{dt} + \frac{R}{L}C e^{-\frac{R}{L}t} = \frac{E}{L}\sin(\omega t)$ (A-2)

$(e^{-\frac{R}{L}t} )\frac{dC}{dt} = \frac{E}{L}\sin(\omega t)$
$\frac{dC}{dt} = e^{\frac{R}{L}t}\frac{E}{L}\sin(\omega t)$

${\rm dC} = e^{\frac{R}{L}t}\frac{E}{L}\sin(\omega t){dt}$

両辺を積分すれば、

$\int {\rm dC} = \int e^{\frac{R}{L}t}\frac{E}{L}\sin(\omega t){dt}$
$\int {\rm dC} = \frac{E}{L}\int e^{\frac{R}{L}t}\sin(\omega t){dt}$

$\int e^{\frac{R}{L}t}\sin(\omega t){dt}$
は部分分数法で展開できて

$\int e^{\frac{R}{L}t}\sin(\omega t){dt}$
$=\int (\frac{L}{R}e^{\frac{R}{L}t})'\sin(\omega t){dt}$
$= \frac{L}{R}e^{\frac{R}{L}t}\sin(\omega t) - \int (\frac{\omega L}{R}e^{\frac{R}{L}t})\cos(\omega t){dt}$
$= \frac{L}{R}e^{\frac{R}{L}t}\sin(\omega t) - \frac{\omega L}{R}\int (e^{\frac{R}{L}t})\cos(\omega t){dt}$

次に
$\int (e^{\frac{R}{L}t})\cos(\omega t){dt}$
$=\int (\frac{L}{R}e^{\frac{R}{L}t})'\cos(\omega t){dt}$
$= \frac{L}{R}e^{\frac{R}{L}t}\cos(\omega t) + \int (\frac{\omega L}{R}e^{\frac{R}{L}t})\sin(\omega t){dt}$

次に
$I=\int e^{\frac{R}{L}t}\sin(\omega t){dt}$ とおくと

$I = \frac{L}{R}e^{\frac{R}{L}t}\sin(\omega t) -\frac{\omega L}{R}\Big(\frac{L}{R}e^{\frac{R}{L}t}\cos(\omega t) + \frac{\omega L}{R} I \Big)$
$I = \frac{L}{R}e^{\frac{R}{L}t}\sin(\omega t) -\frac{\omega L^2}{R^2}e^{\frac{R}{L}t}\cos(\omega t) - \frac{(\omega L)^2}{R^2} I$

$\frac{R^2 + (\omega L)^2}{R^2}I = \frac{L}{R}e^{\frac{R}{L}t}\sin(\omega t) -\frac{\omega L^2}{R^2}e^{\frac{R}{L}t}\cos(\omega t)$
$I = \frac{RL}{R^2 + (\omega L)^2}e^{\frac{R}{L}t}\sin(\omega t) -\frac{\omega L^2}{R^2 + (\omega L)^2}e^{\frac{R}{L}t}\cos(\omega t)$

$I = \frac{L}{R^2 + (\omega L)^2}e^{\frac{R}{L}t}\Big( R\sin(\omega t)-\omega L\cos(\omega t)\Big)$

ゆえに

$\int {\rm dC} = \frac{E}{R^2 + (\omega L)^2}e^{\frac{R}{L}t}\Big( R\sin(\omega t)-\omega L\cos(\omega t)\Big)$

$C = \frac{E}{R^2 + (\omega L)^2}e^{\frac{R}{L}t}\Big( R\sin(\omega t)-\omega L\cos(\omega t)\Big) + D$
(は0を含む任意の実数(積分変数)) (A-5)

■ 変数を元に戻す

(A-2)に(A-5)を代入すると
$u =(\frac{E}{R^2 + (\omega L)^2}e^{\frac{R}{L}t}\Big( R\sin(\omega t)-\omega L\cos(\omega t)\Big) + D) e^{-\frac{R}{L}t}$

$u =\frac{E}{R^2 + (\omega L)^2}\Big( R\sin(\omega t)-\omega L\cos(\omega t)\Big) + De^{-\frac{R}{L}t}$

$u = i^{(1)}$ (3-2)とおいて
u = i なので
$i(t) = \frac{E}{R^2 + (\omega L)^2}\Big( R\sin(\omega t)-\omega L\cos(\omega t)\Big) + De^{-\frac{R}{L}t}$

$= De^{-\frac{R}{L}t} + \frac{E}{\sqrt{R^2 + (\omega L)^2}}\Big( \frac{R}{\sqrt{R^2 + (\omega L)^2}}\sin(\omega t)-\frac{\omega L}{\sqrt{R^2 + (\omega L)^2}}\cos(\omega t)\Big)$

$= De^{-\frac{R}{L}t} + \frac{E}{\sqrt{R^2 + (\omega L)^2}}\Big(\cos(\varphi)\sin(\omega t)-\sin(\varphi)\cos(\omega t)\Big)$
$= De^{-\frac{R}{L}t} + \frac{E}{\sqrt{R^2 + (\omega L)^2}}\Big(\sin(\omega t)\cos(\varphi)-\sin(\varphi)\cos(\omega t)\Big)$

$= De^{-\frac{R}{L}t} + \frac{E}{\sqrt{R^2 + (\omega L)^2}}\sin(\omega t - \varphi)$

ただし、 $\varphi = \arctan(\frac{\omega L}{R})$

ここで、 t=0 のとき i(0)=0 の初期条件の元では、

$0 = D + \frac{E}{\sqrt{R^2 + (\omega L)^2}}\sin(- \varphi)$
$0 = D - \frac{E}{\sqrt{R^2 + (\omega L)^2}}\sin(\varphi)$
$D = \frac{E}{\sqrt{R^2 + (\omega L)^2}}\sin(\varphi)$

$i(t)= \frac{E}{\sqrt{R^2 + (\omega L)^2}}\sin(\varphi)e^{-\frac{R}{L}t} + \frac{E}{\sqrt{R^2 + (\omega L)^2}}\sin(\omega t - \varphi)$

$i(t)= \frac{E}{\sqrt{R^2 + (\omega L)^2}}\Big(\sin(\varphi)e^{-\frac{R}{L}t} + \sin(\omega t - \varphi)\Big)$

ソースコード

$\frac{di(t)}{dt} = \frac{\sin(\omega t + \theta) - Ri(t)}{L}$
の
$f(t,i) = \frac{\sin(\omega t + \theta) - Ri(t)}{L}$
として

$k_{1}(n)=f(t_n , i_n )$ (1-1)
$k_{2}(n)=f(t_n+\frac{h}{2},i_n+\frac{h}{2}*k_{1}(n))$ (1-2)
$k_{3}(n)=f(t_n+\frac{h}{2},i_n+\frac{h}{2}*k_{2}(n))$ (1-3)
$k_{4}(n)=f(t_n+h,i_n + k_{3}(n)*h)$ (1-4)

$i_{n+1} = i_n + \frac{1}{6}( k_{1}(n) + 2k_{2}(n) + 2k_{3}(n) + k_{4}(n) )*h$ (1-5)

// 4ziRungeKuttaSample.cpp : コンソール アプリケーションのエントリ ポイントを定義します。
//

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

#define PI 3.1415926

const static double Em=1.0;                              //電圧源   
const static double R=1.0;                               //抵抗
const static double L=1.0;                               //コイル
const static double freq=0.1;                            //周波数0.1Hz
const static double omega=2*PI*freq;                     //角周波数
const static double theta=-PI/4.0;
const static double phi=atan(omega*L/R);
const static double tau = L/R;                           //時定数
const static double dh=(1/freq)/30.;                     //1周期の30分割
const static double Im=1.0/sqrt(R*R+(omega*L)*(omega*L));

double ii4(double t,double cur);

int _tmain(int argc, _TCHAR* argv[])
{
 
	errno_t err; 
	FILE *fp_output;
	
	//変数
	double t=0;    //実時間
	double e=0; //初期電圧
	double i1=0.0;       //初期電流(厳密解)
	double i4=0.0;       //初期電流(4次のルンゲクッタ法)
	double k[4]={0,0,0,0};
	
	if( err = fopen_s(&fp_output,"data.txt","w")  != 0){ exit(2);}
	
	t=0;
	for(int n=1; n< 3*(int)((1/freq)/dh); n++){
		
		//これは t=n*dh(n=1,2,3,...)の計算となる
		k[0]=ii4(t,i4);
		k[1]=ii4(t+dh/2,i4+dh/2*k[0]);
		k[2]=ii4(t+dh/2,i4+dh/2*k[1]);
		k[3]=ii4(t+dh,i4+dh*k[2]);

		i4=i4+(1./6.)*(k[0]+2*k[1]+2*k[2]+k[3])*dh;
		
		//n=1から始まるから t=n*dh(n=1,2,3...)
		t=(double)n*dh;

		//電源電圧
		e=Em*sin(omega*t+theta);
		//理論値電流
		i1=Im*(sin(omega*t+theta-phi)-exp(-t/tau)*sin(theta-phi));

		printf("t = %f e=%f i1 = %f i4 = %f\n",t,e,i1,i4);
		fprintf_s(fp_output,"%f %f %f %f\n",t,e,i1,i4);
	}

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


double ii4(double t,double cur){

	double calcCur=(sin(omega*t+theta)-R*cur)/L;
	return calcCur;
}

Matlabによる可視化

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5.3 4次のルンゲクッタ法
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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('電圧',...,
'厳密解',...,
'4次のルンゲクッタ法',...,
1);
set(h,'FontSize',15)

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

print -djpeg rungekutta4_method.jpg

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

#define PI 3.1415926

const static double Em=1.0;                              //電圧源   
const static double R=1.0;                               //抵抗
const static double L=1.0;                               //コイル
const static double freq=0.1;                            //周波数0.1Hz
const static double omega=2*PI*freq;                     //角周波数
const static double theta=-PI/4.0;
const static double phi=atan(omega*L/R);
const static double tau = L/R;                           //時定数
const static double dh=(1/freq)/30.;                     //1周期の30分割
const static double Im=1.0/sqrt(R*R+(omega*L)*(omega*L));

double func1(double t,double cur);
double rungeKutta4zi(double t,double i4);

int _tmain(int argc, _TCHAR* argv[])
{
  
	errno_t err;  
	FILE *fp_output;
	
	//変数
	double t=0;          //実時間
	double e=0;          //初期電圧
	double i1=0.0;       //初期電流(厳密解)
 	double i4=0.0;       //初期電流(4次のルンゲクッタ法)
	double i3=0.0;       //初期電流(3次のルンゲクッタ法)
	
	
	if( err = fopen_s(&fp_output,"renritu_bibun_kekka.txt","w")  != 0){ exit(2);}
	
	t=0;
	for(int n=1; n< 3*(int)((1/freq)/dh); n++){
		
	
		i4=rungeKutta4zi(t,i4);
 		i3=rungeKutta4zi(t,i3);
		//n=1から始まるから t=n*dh(n=1,2,3...)
		t=(double)n*dh;

		//電源電圧
 		e=Em*sin(omega*t+theta);
		//理論値電流
 		i1=Im*(sin(omega*t+theta-phi)-exp(-t/tau)*sin(theta-phi));

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

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


double rungeKutta4zi(double t,double i4){
	
	double k[4]={0,0,0,0};
	//これは t=n*dh(n=1,2,3,...)の計算となる
 	k[0]=func1(t,i4);
	k[1]=func1(t+dh/2,i4+dh/2*k[0]);
	k[2]=func1(t+dh/2,i4+dh/2*k[1]);
	k[3]=func1(t+dh,i4+dh*k[2]);

	i4=i4+(1./6.)*(k[0]+2*k[1]+2*k[2]+k[3])*dh;
	
	return i4;
} 

double rungeKutta3zi(double t,double i3){
	
 	double k[4]={0,0,0,0};
 	//これは t=n*dh(n=1,2,3,...)の計算となる
	k[0]=func1(t,i3);
	k[1]=func1(t+dh/2,i3+dh/2*k[0]);
	k[2]=func1(t+dh,(-k[0]+2*k[1])*dh);

	i3=i3+(1./6.)*(k[0]+4*k[1]+k[2])*dh;
	
	return i3;
} 

double func1(double t,double cur){
 
	double calcCur=(sin(omega*t+theta)-R*cur)/L;
	return calcCur;
}

Matlabで描画

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3 4次のルンゲクッタ法
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'renritu_bibun_kekka.txt' 

t=renritu_bibun_kekka(:,1);
v1=renritu_bibun_kekka(:,2);
ir=renritu_bibun_kekka(:,3);
i3=renritu_bibun_kekka(:,4);
i4=renritu_bibun_kekka(:,5);

plot(t,v1,'r-','linewidth',2);
hold on
plot(t,ir,'b-','linewidth',2);
hold on
plot(t,i3,'g+','linewidth',2);
hold on
plot(t,i4,'mo','linewidth',2);
hold on

grid on

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

h=legend('電圧',...,
'厳密解',...,
'3次のルンゲクッタ法',...,
'4次のルンゲクッタ法',...,
1);
set(h,'FontSize',15)

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

print -djpeg rungekutta4and3_method.jpg