日記帳だ! with Tux on Libserver

二度目の大改造!! 日記帳…か?を継承し、より柔軟でパワフルなBlogに変身しました。

RSSに対応しています。リンク・コメント・トラックバックは自由にしていただいてほぼ問題ありません。
RSS購読方法、僕のリンク・コメント・トラックバックについての考えを読むことをおすすめします。

<< 過去

未来 >>

高速で複素積分するC++

数学のノートやらをチェックしていたのだが、どうも授業中に解いた複素関数の積分を間違っていたらしい。
随分へんてこなミスだったので気づけたが、なかなか気づくことは難しい。
普通の積分ならば不定積分に値を投げ込んだら終わりなので、不定積分を微分してチェックするとかできるが、
複素関数の場合ルートによっても変わってくるので楽な話ではない。
複素関数の積分というのはルートにより変わってくるものです。
まぁ変わらない場合もあるそうなのですが変わることは十分にあると。
どうやって積分するのかというと、経路を表す関数を用意します。z(t)などとしておきましょうか。
tっていうのは実数ね。それで f(z(t))って関数にしてtで置換積分するという感じですね。
これでうまいことできます。それさえ気をつければ複素関数だからといって特別なことはないはず。
まぁそりゃ複素数が出てきたりはしますけどね。
そこで数値積分しちゃいましょう。ちょうどC++にはstd::complexってやつもあるんですし。
#include <complex>
#include <iostream>
#include <cmath>
const double PI=3.14159265258979;
typedef std::complex<double> DComplex;
struct CirclePath{
double begin,end,delta;
DComplex operator()(double t){
return DComplex(std::cos(t),std::sin(t));
}
CirclePath(double begin,double end,double delta):begin(begin),end(end),delta(delta){}
};
template <typename TFunc,typename TPath>
DComplex Integral(TFunc func,TPath path){
DComplex result(0,0);
DComplex z1,z2;
bool first=true;
for(double t=path.begin;std::fabs((t-path.end)/path.delta)>1;t+=path.delta){
z2=path(t);
if(!first) result+=func(z1)*(z2-z1);
else first=false;
z1=z2;
}
return result;
}
int main(void){
auto f=[](const DComplex& z)->DComplex{return 1.0/(z*z+9.0);};
std::cout << Integral(f,CirclePath(-PI/2,PI/2,1e-5)) << std::endl;
return 0;
}
いちいちdoubleで構成された複素数型ということでstd::complex<double>と書くのが本来の所だが、
さすがにこれは僕でも嫌気がさしてきたのでtypedefでDComplexとしておきました。
まず経路オブジェクト、さっきのz(t)を定義するものですね。
beginにスタート、endにゴール、deltaに刻み幅を入れた関数オブジェクトにしてあります。
あと複素関数オブジェクトだが、これはC++0xのつもりなのでラムダ式を使うことにした。
Integral関数はまじめに積分しているだけ。簡単ですね。
ここではf(z)=1/(z2+9)という関数を半径1の円の-iのところからiのところに反時計回りに回っている。
これを数学的にまじめに計算しようとすると難しい。
けどコーシーの積分定理を用いれば-iからiまで一直線に積分した結果と等しいということになるので簡単になる。
その結果は i(loge2)/3=0.2310i となる。
けど計算ミスしてて、 (2i/3)tan-1(1/3)=0.2150i になるとか書いてあった。
実は -1から1からにすると結果が (2/3)tan-1(1/3)=0.2150 になるんですが、この辺と混乱してしまったらしい。
ところがすぐに間違いに気づけなかったというマヌケな話。
けどこのプログラムを使うと答えがすぐにわかる。
(1.67633e-006,0.231049)など出てくる。実部は非常に小さいので無視、結果は0.2310iですね。
確かに合っていることがわかる。
コンピュータというのは本当にありがたい。
数字でしか出せませんがなかなか良い結果が出てきますね。
しかもC++で作ったプログラムは驚くほど高速で動きます。
コンパイルの時間を入れてもすぐですね。なので非常にお得。
初めこれをPerlかRubyで書こうと思ったんだが時間かかるだろうなと思ってやめたんだよね。
これである程度の答え合わせはできますね。
Author : Hidemaro
Date : 2009/09/20(Sun) 23:54
電気・数学・物理 | Comment | trackback (0)
blog comments powered by Disqus

トラックバック

トラックバックURL取得

Tools