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

参考文献及び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) 微分方程式の解を数値計算で求める

(c) 積分の説明

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

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

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

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

はじめに

常微分方程式を解く前に、近似とテーラー展開を復習する

近似

今、 x = a のときの関数値 f(a) が求められるとするとき、
変数が a からわずかに違う x=a+\Delta xになったときの
f(a+\Delta x)は？

f(a+\Delta x) \approx f(a)+ (\frac{df(a)}{dx})\Delta x (1-1)

[重要]

f(a) が分かっているとき
f(a)と f(a+\Delta x)を結ぶ直線を
aにおける曲線の接線と等しいとみなす。

(例題) 1辺10cmの立方体が熱されて1辺が1mm伸びたときその体積は？

[回答]

一辺をxとし、体積をyとすると

y=x^3

10^3 = 1000 \rm{cm}^3

a=10とし、aから \Delta x=0.1離れた点を考える。

○近似式で解く

\frac{dy(10)}{dx} = 3x^2 =3 \times 10^2 = 300

f(a+\Delta x) \approx f(a)+ (\frac{df(a)}{dx})\Delta x (1-1)
y(10+0.1) = f(10) + \frac{dy(10)}{dx} \times 0.1
=1000 + 300 \times 0.1 = 1030 \rm{cm}^3

○正確な値は

10.1^3 = 1030.301 \rm{cm}^3

○誤差率は

\sigma = \frac{ 1030.301 - 1030 }{1030.301} = 0.029%

テーラー展開

ある点aにおける関数f(a)が決まるとき、
ある曲線は、近似を1次(直線),2次(曲線),3次...というように
無限個の近似曲線の集合によって、もとの曲線に近づけることができる。
これをテーラー展開という。

f(a+x) = f(a) + \frac{f'(a)}{1!}x + \frac{f''(a)}{2!}x^2 + \frac{f'''(a)}{3!}x^3 + \cdots

= f(a) +\sum^{\infty}_{n=1}\frac{f^{n}(a)x^{n}}{n!}

○例

e^x = 1+x+\frac{x^2}{2!} + \frac{x^3}{3!} + \cdots
\log(1+x) = x - \frac{x^2}{2!} + \frac{x^3}{3!}- \frac{x^4}{4!} +\cdots

テイラー法

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

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

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

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

これを初期条件

y(x_0) = y_0

いま、 x_0 の点からn \cdot h離れた点
( n = 1,2,\cdots r+1は整数で、 nh < 1とする)

の曲線の方程式を求める。

y_n = y(x_0 + nh)
= y_0 + \frac{y_{0}'}{1!}\cdot(nh) + \frac{y_{0}''}{2!}\cdot(nh)^2 + \frac{y_{0}'''}{3!}\cdot(nh)^{3} + \cdots + \frac{y_{0}^(r)}{r!}\cdot(nh)^r + \frac{y_{0}^(r+1)}{(r+1)!}\cdot(nh)^{(r+1)} (2-1)

ここで、

y' = \frac{dy}{dx} = f(x,y(x))

なので

y'' = \frac{d^2y}{dx^2}
= \frac{3$\partial~f(x,y)}{3$\partial~x} + \frac{3$\partial~f(x,y)}{3$\partial~y}\frac{3$dy}{3$dx}
= f_x(x,y) + f_y(x,y)y'

同じように
y'''= \frac{d}{dx}\Big(f_x(x,y) + f_y(x,y)y'\Big)
= \frac{d}{dx}f_x(x,y) + \frac{d}{dx}\Big(f_y(x,y)y'\Big) (2-2)

(2-2)の第一項は、

\frac{d}{dx}f_x(x,y)
= f_{xx}(x,y) + \frac{\partial~f_{x}(x,y)}{\partial~y}\frac{dy}{dx} = f_{xx}(x,y) + f_{xy}(x,y)y'

(2-2)の第二項は、
\frac{d}{dx}\Big(f_y(x,y)y'\Big)
= \frac{d}{dx}\Big(f_y(x,y)\Big)y' + \frac{d}{dx}(y')\Big(f_y(x,y)\Big) (2-3)

(2-3)の第一項は
\frac{d}{dx}\Big(f_y(x,y)\Big)y'
= \Big( f_{xy}(x,y) + \frac{\partial~f_{y}(x,y)}{\partial~y}\frac{dy}{dx} \Big)y'
= \Big( f_{xy}(x,y) + f_{yy}(x,y)y' \Big)y'
= f_{xy}(x,y)y' + f_{yy}(x,y)y'^2

(2-3)の第二項は

\frac{d}{dx}(y')\Big(f_y(x,y)\Big)
= y'' f_y(x,y) + y' \frac{d}{dx}\Big(f_y(x,y)\Big)
= y'' f_y(x,y) + y' \Big(f_{xy}(x,y) + \frac{\partial~f_{y}(x,y)}{\partial~y}\frac{dy}{dx}\Big)
= y'' f_y(x,y) + y' \Big(f_{xy}(x,y) + f_{yy}(x,y)y'\Big)
= y'' f_y(x,y) + y'f_{xy}(x,y) + y'^{2}f_{yy}(x,y)

ゆえに

y'''= \frac{d}{dx}\Big(f_x(x,y) + f_y(x,y)y'\Big)
= \frac{d}{dx}f_x(x,y) + \frac{d}{dx}\Big(f_y(x,y)y'\Big)
= f_{xx}(x,y) + f_{xy}(x,y)y' + f_{xy}(x,y)y' + f_{yy}(x,y)y'^2 + y'' f_y(x,y) + y'f_{xy}(x,y) + y'^{2}f_{yy}(x,y)

= f_{xx}(x,y) + 3f_{xy}(x,y)y' + 2f_{yy}(x,y)y'^2 + y'' f_y(x,y)
([再考] この計算あってるかな？？)

このようにy',y'',y''',\cdotsを順次計算し、初期条件x_0のときの結果を
(2-1)に代入するとx_0の点からx_0 < x < x_0+1の範囲の yを求めることができる。

xがx_0 + 1 < xになってしまったら、
超える直前の xを x_1とし
再度1を超えないように計算をすればよい。

[例題1] \frac{dy}{dx} = -\sin3x -3y, y(0)=y_0=1 をテイラー法でとく

[条件]
x \lt 1 ( nh \lt 1)
h = 1/40
厳密解 y = \frac{1}{6}(\cos(3x) - \sin(3x) + 5e^{-3x}

y' = -\sin3x - 3y
y'' = -3\cos3x - 3y'
y''' = 9\sin3x - 3y''
y^{(4)} = 27\cos3x - 3y^{(3)}
y^{(5)} = -81\cos3x - 3y^{(4)}

\cdots

初期値 y(0)=y_0=1 を代入すると

y(0) = 1
y'(0) = -\sin3x - 3y = 0 -3= -3
y''(0) = -3\cos3x - 3y' = -3 +9 = 6
y'''(0) = 9\sin3x - 3y'' = -18
y^{(4)}(0) = 27\cos3x - 3y^{(3)} = 27 + 54 = 81
y^{(5)}(0) = -81\sin3x - 3y^{(4)} = -243

\cdots

いま、 h = 1/40 とすると

y_n = y(x_0 + nh)
= y_0 + \frac{y_{0}'}{1!}\cdot(nh) + \frac{y_{0}''}{2!}\cdot(nh)^2 + \frac{y_{0}'''}{3!}\cdot(nh)^{3} + \cdots + \frac{y_{0}^(r)}{r!}\cdot(nh)^r + \frac{y_{0}^(r+1)}{(r+1)!}\cdot(nh)^{(r+1)} (2-1)

= 1 + \frac{-3}{1}\cdot(\frac{n}{40}) + \frac{6}{2}(\frac{n}{40})^2 + \frac{-18}{6}(\frac{n}{40})^{3} + \frac{81}{24}(\frac{n}{40})^{4} + \frac{-243}{120}(\frac{n}{40})^{5} + \cdots

解は5項までもとめた。テイラー法の項数を増やせばより厳密解に近くなる。

以下のプログラムの結果、テイラー法の項数は10ぐらいあれば厳密解とほぼ等しくなる

[例題2] [例題1]をプログラム化する

#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 n   = 0;
	double h   =1/40.;
	double x0  =0;
	double y2=0,y5=0,y10=0,x=0,correcty=0;
	double yy[10];

	yy[0]= 1;                          //0階微分 x=0のときのy(x)
	yy[1]= -sin(3*x) - 3*yy[0];        //1階微分 x=0のときのy'(x)
	yy[2]= -3*cos(3*x)-3*yy[1];        //2階微分 x=0のときのy''(x)
	yy[3]=  pow(3.,2.)*sin(3*x)-3*yy[2];
	yy[4]=  pow(3.,3.)*cos(3*x)-3*yy[3];	
	yy[5]= -pow(3.,4.)*sin(3*x)-3*yy[4];	
	yy[6]= -pow(3.,5.)*cos(3*x)-3*yy[5];
	yy[7]=  pow(3.,6.)*sin(3*x)-3*yy[6];	
	yy[8]=  pow(3.,7.)*cos(3*x)-3*yy[7];
	yy[9]= -pow(3.,8.)*sin(3*x)-3*yy[8];


	if( err = fopen_s(&fp_output,"data.txt","w")  != 0){ exit(2);}
	
	for(float n=1; (n*h) < 1;n++){
		x=x0+n*h;

		y2 = yy[0] 
			+ (1./kai(1))*yy[1]*pow((n*h),1) 
			+ (1./kai(2))*yy[2]*pow((n*h),2);

		y5 = yy[0] 
			+ (1./kai(1))*yy[1]*pow((n*h),1) 
			+ (1./kai(2))*yy[2]*pow((n*h),2)
			+ (1./kai(3))*yy[3]*pow((n*h),3)
			+ (1./kai(4))*yy[4]*pow((n*h),4)
			+ (1./kai(5))*yy[5]*pow((n*h),5);

		y10 = yy[0] 
			+ (1./kai(1))*yy[1]*pow((n*h),1)  
			+ (1./kai(2))*yy[2]*pow((n*h),2) 
			+ (1./kai(3))*yy[3]*pow((n*h),3) 
			+ (1./kai(4))*yy[4]*pow((n*h),4) 
			+ (1./kai(5))*yy[5]*pow((n*h),5) 
			+ (1./kai(6))*yy[6]*pow((n*h),6) 
			+ (1./kai(7))*yy[7]*pow((n*h),7) 
			+ (1./kai(8))*yy[8]*pow((n*h),8) 
			+ (1./kai(9))*yy[9]*pow((n*h),9) ;

		correcty=f(x);

		printf("x = %f cy = %f y2 = %f y5 = %f y10 = %f\n",x,correcty,y2,y5,y10);
		fprintf_s(fp_output,"%f %f %f %f %f\n",x,correcty,y2,y5,y10);
	}

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

double kai(int in){

	double ans=1;
	for(int n=in;n > 0;n--){
		ans=ans*(double)n;
	}
	return ans;
} 

double f(double x){
	double y;
	y=(1./6.) * (cos(3*x)-sin(3*x)+5*exp(-3*x));
	return y;
}

[出力] Matlabスクリプト

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5.1 テイラー法
% written by embedded.samurai
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
echo off
clear all
close all

load 'data.txt'

x=data(:,1);
cy=data(:,2);
y2=data(:,3);
y5=data(:,4);
y10=data(:,5);

plot(x,cy,'r-','linewidth',6);
hold on
plot(x,y2,'b+-','linewidth',2);
hold on
plot(x,y5,'g+-','linewidth',2);
hold on
plot(x,y10,'m+-','linewidth',2);
hold on
grid on

xlabel('x','Fontsize',18)
ylabel('y','Fontsize',18)

h=legend('厳密解',...,
'テイラー法(第2項まで)',...,
'テイラー法(第5項まで)',...,
'テイラー法(第10項まで)',...,
3);
set(h,'FontSize',15)

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

print -djpeg telor_method.jpg
Generated by ひとりWiki version 1.15.1