ArchLinux,MATLAB,制御理論など(の予定)

数値での微分についてちょっと考える

2017-02-11

今更だけど数値計算で微分を表現することについて考えてみる。

動機 #

最近数値的に微分計算をしたいと思うことがあるが、 一体誤差がどれぐらいで見積もれるのか理解しておくべきだと思ったのでまとめておく。
今更感が丸出しだけど。。。

問題設定 #

時間間隔\(dt\)でデータ列\(f_k(k=0,1,...,K)\)が得られたとする。この時、それぞれの時刻での微分値が知りたい。

線形近似 #

連続→離散 #

当たり前だが連続値の微分はコンピュータ上では表現できない。
そこで連続的な微分を離散的に扱わなければならないが、表現の仕方によって誤差が生じる。

基本はテイラー展開 #

ある連続関数\(f(t)\)が存在していたとき、ある点\(f(t+dt)\)は以下のようにテイラー展開することが出来る(実際にはテイラー級数が収束する場合にのみテイラー展開を定義できるがとりあえず収束するとする)

\[ f(t+dt)=\sum_{n=0}^{\infty}\frac{f^{n}(t)}{n!}dt^n \]

なぜこの形にテイラー展開出来るのか、またテイラー展開の承継などはネットにゴロゴロ転がっているので良いとして、ここからどうやって離散にしていくのか。

\(f(kdt)\)\(f_k\) #

まぁ順当に考えて時間間隔で離散化されているわけだから、離散化する変数を\(k\)として テイラー展開の式をちょっと書き直してみる。
\(t=kdt\)として見ると以下のように書き直せる。

\[ f_{k+1}=\sum_{n=0}^{\infty}\frac{f^{n}_k}{n!}dt^n \]

線形近似なので微小量\(dt\)が2乗以上で効いてくる項をぶち消すと、

\[ f_{k+1}=f_{k}+f'_kdt \]

ここで知りたかったのは微分値なので微分値について解くと

\[ f'_k=\frac{f_{k+1}-f_{k}}{dt} \]

まぁそうだよね!ありきたりな結果!知ってた!

結局この式が何を示しているかというと、ある点に対して次の点との増分を時間間隔 割っているわけなのでまさしく線形的な傾きを求めていることになる。

しかし、線形で近似するともちろん誤差がでかくなる。 そこで、2次曲線で近似してみたらどうだろうというのが次のコンセプト。

2次近似 #

至って簡単なコンセプトで、注目した点\(t_k,f_k\)に対して前後の\(t_{k+1},f_{k+1},t_{k-1},f_{k-1}\)を通るような二次曲線を考え、\(t_k\)における微分値を採択するという方法である。 三点を通る二次曲線を\(g(t)=at^2+bt+c\)とすると、三点の値を用いて以下の式が成り立つ。

\[ A = \left( \begin{array}{ccc} t^2_{k-1} & t_{k-1} & 1 \\ t^2_k & t_k & 1 \\ t^2_{k+1} & t_{k+1} & 1 \end{array} \right)\\ \theta = \left( \begin{array}{c} a \\ b \\ c \end{array} \right)\\ F = \left( \begin{array}{c} f_{k-1} \\ f_k \\ f_{k+1} \end{array} \right)\\ A\theta=F \]

この方程式を解いた結果の\(a,b,c\)\(g'(t)=2at+b\)から

\[ f'_k=2at_k+b \]

と計算すると微分値が得られる。より高次の多項式で近似する場合も同じようにすればいい。

MATLABによる検証 #

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
function dy = diff2(t,y)
n=length(t);
dy=zeros(1,n);
y2=y.*y;
t2=t.*t;
for i=2:n-1
A=det([t2(i-1:i+1)',t(i-1:i+1)',ones(3,1)]);
a=det([y(i-1:i+1)',t(i-1:i+1)',ones(3,1)])/A;
b=det([t2(i-1:i+1)',y(i-1:i+1)',ones(3,1)])/A;
%c=det([t2(i-1:i+1)',t(i-1:i+1)',y(i-1:i+1)'])/A;
dy(i)=2*a*t(i)+b;
end
dy(1)=dy(2);
dy(n)=dy(n-1);
end
dt=0.001;
t=0:dt:1;
y=sin(t)+sin(10*t);
dy=cos(t)+10*cos(10*t);
y1=[y(1),y(1:length(y)-1)];
y2=[y(2:length(y)),y(length(y))];
dy1=(y-y1)/dt;
dy2=(y2-y)/dt;
dy3=(dy1+dy2)/2;
dy4=diff2(t,y);
plot(t,y,t,dy,t,dy1,t,dy2,t,dy3,t,dy4)
e1=sum(abs(dy1-dy))/length(dy)
e2=sum(abs(dy2-dy))/length(dy)
e3=sum(abs(dy3-dy))/length(dy)
e4=sum(abs(dy4-dy))/length(dy)

dy1は後退差分近似、dy2は前進差分近似、dy3は中心差分近似、dy4は2次近似である。
この結果は

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
e1 =
0.0418
e2 =
0.0387
e3 =
0.0095
e4 =
1.6241e-04

といった具合。e4がやっぱり小さい。2次近似の効力って感じ。

まとめ #

じゃぁ、何でもかんでも2次近似を用いればいいかというとそうでもない。例えば、計算コストという点で考えると明らかに2次近似のほうが辛い。今回は行列式を解かせる形だったので余計にコストがかかっているが、仮にこれを紙の上で解いて代入するだけの形にしていたとしても、明らかに線形近似よりも計算量は多いだろう。
また、使っている変数の数(注目する時刻に対して用いるデータの数)はどうだろうか。
前進差分近似、後退差分近似では前後とその時刻の2つ、中心差分近似と2次近似では3つのデータを使う必要がある。
今回は全ての時刻での\(f_k\)が出揃った状態で計算するため考えなくても良かったが、時間発展的に計算したい場合利用できるデータは現在と過去のデータのみであり、今回のプログラムで用いた方法の中で後退差分近似以外は使えない。
結局、計算コストと精度の兼ね合いである。データの雰囲気がある程度予測できる時、適切な次数で近似するというセンスが求められる。
私も勉強中。。。

そのうち時間発展的に計算する方法(微分方程式を解く方法)も書いていきたい。


Blog comments powered by Disqus