数値積分や補間，代数方程式の求解は, 変数ではなく関数を与える必要があります．「関数を与える」ってところが独特ですが，みんな同じなので，ここで覚えておきましょう．

数値積分

関数を与えて, 定積分を計算する数値積分にはいろいろな方法があります．ここにはいくつか書いてありますが，まあ例として二重指数型積分公式を使ってみましょう.

問題例1

とりあえず簡単な定積分を行ってみる:

e^{- U^{2}} + \frac{1}{2 \sqrt{π}} \int_{0}^{2 π} \exp [- U^{2} \sin^{2} θ] erf (U \cos θ) U \cos θ d θ = 1

なんてどうでしょう. なお

erf (x) = \frac{2}{\sqrt{π}} \int_{0}^{x} e^{- t^{2}} d t

でしたね．こりゃ楽すぎるかも

問題例2

とりあえず面倒な積分も行ってみる:

\frac{2}{π} \int_{0}^{1 / 2} \frac{1}{\sqrt{x (1 - x)}} d x = 1

とかどうですかね．

オイラは自作派

自作する手もあります．台形則なんて素晴らしいですよね. ええっとこれは，確か台形の面積が上辺＋下辺かける高さ割る2，だったので，区間 $[a b]$ を $N$ 等分して

\int_{a}^{b} f (x) d x = \sum_{i = 0}^{N - 1} \frac{f (x_{i}) + f (x_{i + 1})}{2} Δ x

ってわけですよね．

figs

問題例1

いくらなんでもちょろいですね．まあこんな感じか. まず被積分関数を定義して

double integrand(double th,double U)
{
  auto UC=U*std::cos(th);// auto でもdoubleでも大差ないな
  auto US=U*std::sin(th);// std:: は省略しても同じ
  auto ERF=std::erf(UC); // 誤差関数は標準関数です
  auto EXP=std::exp(-US*US);
  auto COEF=M_2_SQRTPI/4.;  //【Jump to Definition】すれば定義がわかる
  return COEF*EXP*ERF*UC;
}

んで，足し合わせるだけですね：

double sum=0.;
for(int i=0;i<N;i++)
{
  double y0=integrand(dth*i,U);  //細かく言えばy0は前回のy1なので手抜きできる
  double y1=integrand(dth*(i+1),U);
  sum+=dth/2*(y0+y1);
}
sum+=std::exp(-U*U)-1.;

結果を見てみると, $U = 1$ ならばこんな感じ:

$N$	10	20	40	80
誤差	$4.2 \times 10^{- 5}$	$1.9 \times 10^{- 12}$	$- 2.2 \times 10^{- 16}$	$- 1.1 \times 10^{- 16}$

最初は調子良く誤差が減るのに, 最後はあまり変わりません．double 浮動小数点数は，およそ15桁しか有効数字がないので， $10^{- 12}$ レベルから先は計算が怪しげになるからです.

台形則でも十分ですが，シンプソン公式を用いると，もっと良いです. プログラムしてみましょう．

問題例2

では問題例2はどうだろう? こいつは $x \to 0$ で $f (x) \to \infty$ なので積分しにくいですね．変数変換したらもっと簡単になるのはもちろんですが，実戦では，「実は被積分関数は数値データなので，式が不明である．端点で発散することは知っている」「被積分関数の1点を計算するのに千円かかる．なお予算は5万円である」「変数変換はグルに禁止されている」ということが多々発生するわけです．

仕方ないですね． $x \in [ϵ, 1 / 2]$ の積分で誤魔化すことにしましょう. やってみると，まあ誤差はこんな感じですね．

	$N = 10^{4}$	$N = 10^{6}$	$N = 10^{8}$	$N = 10^{10}$
$ϵ = 10^{- 8}$	$0.15$	$0.00092$	$- 0.00013$	$- 0.00013$
$ϵ = 10^{- 10}$	$1.6$	$0.015$	$0.000092$	$- 0.0000010$

天文学的な格子点数が必要です．やはり $x \sim 0$ で等分割ではダメでしょうね．．．といって工夫し始めると，どんどん時間がかかっていきます．

オイラは他力本願

めんどくさくなってきたので，プロが書いたライブラリーに積分してもらいましょう．ライブラリーには，変数ではなくて関数を渡さなければならないので，ちょいと覚えるものがあります．

ラムダ式

ラムダ式とは, 名前がなく, 名前で呼び出すことが不可能な関数です．名前はありませんが，「あれやって，これやって，・・・」という手順はあるわけです．そこで，それを変数に入れると，呼び出せるようになります．なにぶん「手順」という値ですので，

他の関数の引数に入れられる
コピーしたり，関数の値として呼び出し元に返したり

まるで数値のように関数が使えるわけです．ラムダ式の定義は簡単で

auto F=[](){std::out<<"はろー" << std::endl;};

最初の[]がラムダ関数の開始命令, 関数の引数がこれ, return文がなければvoid型ってわけです．値を返却する場合は, return文を入れると, まあ推定してくれます:

auto F=[](int j){std::out<<"はろー" << std::endl;return j+1;};

この例ですと, jが整数 int ってこた j+1も普通 int なので, たぶん，Fの値は int になります．「たぶん」が嫌なら,

auto F=[](int j)->auto{std::out<<"はろー" << std::endl;return j+1;};

あれ？これでは auto 推定なので多分のままだわ. じゃ，

auto F=[](int j)->int{std::out<<"はろー" << std::endl;return j+1;};

これなら，疑いの余地なく int ですね! 最初の[]に微妙に疑問を覚える人は，ここをみてね♡

これらを受け取った変数Fは，数値ではなく，手順を覚えています．どう使うのかといえば,

auto ret=F(); //retには戻り値のコピーが代入される

これで呼び出すことができます．FをPoiPoiという変数にコピーしたら，同じ関数の呼び出しが

auto R2=PoiPoi() ;

になります．なんといっても変数ですから，他の関数の引数に与えることができます.

some_function(other values,....., PoiPoi,....);
あるいは
some_function(other values,....., [](){std::out<<"はろー" << std::endl;return 0;},...);

1行目はまあわかりますね．2行目は，もう面倒くさいので引数の値としてラムダ式を書いています．

では，ラムダ式を使って, ここで追加したライブラリー intDE で積分してみましょう.

問題例2

#include "intde.hpp"
...
class Pakapaka   //被積分関数には，多くの場合パラメータが必要です．それを入れておくクラス
{                //まあ好きなように定義してね
public:
  int count=0;
};
...
int main()
{
...
Pakapaka parm;  //パラメータが入っているインスタンスが parm
...
IntDE myINT;              //積分屋のインスタンスがmyINT
myINT.Parameter=&parm;    //パラメータのポインターを指定
myINT.Initialize(err_request);    //希望の積分精度を指定する. 1e-8とか
auto v=myINT.Integrate([](double x,void *p) //IntDEの被積分関数の引数は必ずコレ!
  {
    auto& P=*(Pakapaka *)p; //あなたが与えたPakapakaクラス
    P.count++;              //パラメータと言っても，値を変えることも可能
    return 2./sqrt(x*(1.-x))/M_PI;  // 積分関数値をreturnしてね
  },xmin,xmax,err);         //下限, 上限を与え, 推定した精度を取得する
...
}

緑の枠とオレンジの枠に何が書いてあるのかわかりませんが，これらはあとで説明しますね．

この例では, 被積分関数値を何回計算したか知りたかったので, 計算するたびにcount++しています. では, 結果です：

要求精度	$10^{- 4}$	$10^{- 8}$	$10^{- 12}$
推定精度	$2.0 \times 10^{- 4}$	$2.0 \times 10^{- 8}$	$2.0 \times 10^{- 12}$
本当の精度	$4.9 \times 10^{- 7}$	$1.1 \times 10^{- 10}$	$1.0 \times 10^{- 15}$
関数計算回数	$N = 15$	$N = 31$	$N = 50$

プロの仕事とは, このようなものです. 到底かないません.

与えられた関数の数値積分を自分でコーディングする人は，事情がない限り，素人と呼ばれます．

ポインター

緑の枠は，ポインターや参照と呼ばれる文法と関連してます．

ポインターとは，変数の値を記録している場所のことです．コンピュータでは，変数の値をメモリー

ddr4

に保存しています．これは，変数を記録しているところが一列に並んで，アドレスで区別しています．

addr

まあこんな感じですかね．で, プログラム中に count と書けば, 1 とか44とか，値がわかるわけですが，&count と書くと，アドレス「東京都葛飾区若松町東一丁目38-8」がわかるわけです．

一体何の役に立つのか？上の例では，ユーザーが利用するパラメータを

myINT.Parameter=&parm;

で与えています．なぜこんなことをするのかと言うと，

myINTを書いた人が，あんたが今からいくつのパラメータを使うつもりなのか，知るわけがない

ためです．myINTに与えるのは，「ここから先，どこまでかは知らんが，パラメータらしいですぜ」と言う情報だけなので，アドレスを渡しているのです．

ad2

ところで，アドレスも「東京都葛飾区若松町東一丁目38-6」という値を持っていますので，これをコンピュータに記録することもあります．例えばdouble値を持つ変数のアドレスであれば，これは

double*

と書きます．例えば double X 変数があったら,

double* X_point = & X;

で, X_pointにアドレスの値が代入されます．C++では空白は意味がないので, double* と書く人や, double *と書く人がいます．X_point はdouble変数が入っている方向を指差していることになりますが，では，その指の先に書いてある値は,

*X_point

で取得できます・・・この辺で，ほとんどの人の頭はこんがらがってきます．では，以下のプログラムは正しいでしょうか．

double X=3.0, Y=2.0;
double* X_point=&X;
Y=X_point;

だめですね．Yはdouble数値なのですが, X_pointは「数値がある方向」です．「方向」はdoubleの数値ではないので，Yには代入できません．

double X=3.0, Y=2.0;
double* X_point=&X;
Y=*X_point;

であれば, X_point方向に書いてある値はdoubleの値3.0ですので, Yの値が3.0に変わります．

では, サンプルの被積分関数の引数にある

void * p

はなにか？これは，「いやー何が入っているかわからないんだけど，アドレスであることは確か」という表現です．IntDEを書いた人は，ここはユーザーがパラメータを入れるところ(myINT.Parameterで与えたやつ)ということ以外は，情報が分かりませんからね. これ以上は何もかけないわけですよ．IntDE作者の森正武先生はすでに御他界されているので，今から君がなにをするといっても何もできないわけです．

でも，IntDEを使うプログラムを書くユーザーは, p に何のアドレスが入っているかは, わかるわけです．だってほれ，

myINT.Parameter=&parm;

で与えたんですからね．こいつには, 俺様のパラメータクラス Pakapaka のインスタンスparmのアドレスが入っているわけです．で，それを利用するのであれば，「いやー何が入っているかわからないんだけど，アドレスであることは確か」というものから，「それが何かはともかくPakapaka のインスタンスのアドレス」に変換しなければなりません．この作業をキャストと呼びます．サンプルにあるように,

(Pakapaka *)p

と書くと, これはそれが何かはともかくPakapaka のインスタンスのアドレスになります．myINT.Parameter=&parm と設定したことを思い出すと, それは実はparmのアドレスである．というわけで，

auto P=*(Pakapaka *)p;

と書くと, Pはparmのアドレスに書いてある値，つまりにparmの値になる. おっと？かすかに違いますね．サンプルには

auto& P=*(Pakapaka *)p;

と書いてありますよ？

参照

上のサンプルでは，「IntDE様が，関数を何回お呼び出しなされたのか，調べてみたい」ということで，途中でパラメータを変更しています：

P.count++; //パラメータと言っても，値を変えることも可能

なるほど，関数を呼び出すたびに，parmに保存されている count の値を1増やす作戦ですね．

作戦は了承しましたが，これは

auto P=*(Pakapaka *)p;

と書くと, 動作しません．順を追って考えてみましょう．

はじめ，p はポインターであった．
(Pakapaka *)p によって, Pakapaka インスタンス(parm)の方向が得られた．
*(Pakapaka *)p によって, parmに書いてある値となった．その値を，そうですね，VAL と書くことにしましょう．
auto P=VAL を実行した

最後が問題ですね．実は．プログラムの = は等号ではなく，コピー！なのです．だから，PはオリジナルVALのコピーです．で，衝撃的な事実を指摘しておきましょう：「コピーとオリジナルは，実は，別の物体である」

さて，上に注意しつつ，何をしているのか考えてみましょう，

P.count を1増やす

で，プログラム終了直前に，オリジナルであるVALのcountを調べると，呼び出し回数ゲット〜!

・・・おかしいですね．上の手順をよくみると，

VALのコピーをとる
コピーの方に落書きをする
VALに落書きが書いてあるはず！

と期待していることになります．いや〜，これでVALに落書きが書いてあったら，オカルトですよね・・・

これは初心者が陥りやすく，簡単に数週間無駄にしてしまうものです．次の鉄則を覚えましょう．

= はコピーを作成する!
関数で return してもコピーになる!
コピーを書き換えても，本物には影響がない！
だから，コピーするな！

おいおい，どうやってコピーを避けるのか？それを実現するのが参照です．サンプルの

auto& P=*(Pakapaka *)p;

この紫部分が参照です．例えば, double X の参照を定義するのは

double & X_ref

です．ですが参照はコピーではない ＝ を実現しようとするものですから，double & X_ref だけでは機能しません．等号の右辺がないからです．参照を定義するには，必ず値が必要です．つまり

double X=3.0;
double & X_ref=X;

で機能します．X_ref はXを参照しています．よって

X_ref=5.0;

とすると, X_refの値つまりXの値が変更されます．サンプルにある

auto& P=*(Pakapaka *)p;

では, Pはユーザーが myINT.Parameter=&parm で与えたパラメータインスタンスの本物の参照になります．コピーではないです！　だから, P.count を1増やすと, 本物の count が増大するので，計算終了後には, countに呼び出し回数が保存されています．

実はポインターの先にあるクラスのメンバーをゲットすることができます. 上の例では

auto P=(Pakapaka *)p;    // この場合，PはPakapakaインスタンスの方向（ポインター）
P->count++;              // ポインターの先の値のcountメンバーが P->count

これでも無事に動作します．ー＞という記号が，なんとも，うざいけど．

ポインターはC言語で1970年ごろ発明され，強力便利で広がったのですが，とにかくハッキングに弱いという欠点があり，現在では駆逐対象の文法です. C++は古い言語なので，結構利用する局面があるでしょう．

ポインターを利用した多数のPCウイルスに怒り狂ったMicrosoftは, C言語からポインターを取り除き，存在しないものは使えないというC#言語を発案し，それでWindowsを記述できるように.NETを開発しました．Windowsが兎にも角にも動作しているところを見ると，ポインターって，本当に必要ないんですね．でもDelegateがむやみに難しいんだけど

参照は，ポインターの代わりとして，新しいプログラミング言語では広く使われています．

参照が引数の関数

参照を利用すると，関数の引数に与えた変数を変更できます．普通の関数では

int my_func(double x,double y){....}

となっており, xとyは「引数」であって，その値がmy_func()の実行で変化しません．my_func()の実行中に

と行っても，my_func()が保有する「コピー」の値が変化するだけで，呼び出し元のオリジナルxとyの値は変化しません．ところが！

int your_func(double& x, double &y){....}

のように, xとyの参照を引き渡すと話が違います． your_func()の実行中に

x=x+1;

のように，xの値を変化させると・・・なんせ参照なので，呼び出し元のオリジナルxとyの値も変化してしまいます．なにせオリジナルを渡してしまったので，z=your_func(x,y)を実行すると, zが得られるだけではなく，x, yの値が変化するという挙動をしてしまいます．便利ではありますが，

z=your_func(1.0, 2.3);

みたいに，変数ではないもののオリジナルを渡そうとすると「いやそれ変数じゃねえし」とエラーしますので，注意しておいてください．上のIntegrateの引数のerrは，この仕組みを使って，推定精度を呼び出し元に返却しています．

コメントを投稿するにはログインしてください