JacobiのLaplaceは解く！

まずmain()はこんな感じ

#include "myArray.hpp"
#include "jacobi.hpp"
#include <iomanip>
int main(int argc, const char * argv[]) {
  jacobi Solver(11); //ソルバーを初期化. 格子点数10x10
  double e=1.0;      //前回とのずれ. とりあえずバカでっかくしておかないと1回目が動かねえ
  int n=0;           //イタレーション数
  while ((e=Solver.step())>1e-15) //1ステップやって, ズレが\(10^{-15}\)以上ならもう一回だ
  {
    if (++n%100==0)  //イタレーション数1増やして100で割って0なら画面出力
      std::cout << "Step:Diff= " << n << " " << e << std::endl;
  }
  std::cout << "Finish Step:Diff= " << n << " " << e << std::endl;
  return 0;
}

となると, jacobiクラスへの要求は

jacobiクラスは, 格子点数を与えると配列データを作成して，多分初期値とか境界条件とか設定する
jacobi.step()は, 1回進めて, ずれ量の最大値を返却する

必要があるわけですね．では実装しましょう．

#ifndef jacobi_hpp
#define jacobi_hpp
#include "myArray.hpp"
class jacobi
{
  myArray<double> U,U_old;
public:
  jacobi(int N):U(N,N),U_old(N,N)
  {
    for(auto& uv:U) uv=0.0;//初期条件. 境界条件はx=0でU=1, 残りは0ってことで
    for(int i=0;i<U.size(0);i++) U(i,0)=1.0;
  }
  double step()
  {
    double diff=0.;
    U_old=U;
    for(int i=1;i<U.size(0)-1;i++) for(int j=1;j<U.size(1)-1;j++)
    {
      U(i,j)=(U_old(i-1,j)+U_old(i+1,j)+U_old(i,j-1)+U_old(i,j+1))/4.;
      diff=std::max(diff,std::abs(U(i,j)-U_old(i,j)));
    }
    return diff;
  }
};
#endif /* jacobi_hpp */

これくらいですかね.

ここでは自作原理主義配列クラスを使っています．
- ところが問題があります．このmyArray<double>クラスを初期化するには, 配列サイズが必要なのです！ところが「myArray<double> U(ここらへん，未確定なんだけど)」という問題が出現します．
- それには, ここで学んだメンバー変数初期化リストを使えば良いのです！このへんですね．
ここでは, 作成した配列の全要素に，ここで学んだrange based for で初期値を設定しています．
- 実は, 自作クラスには，デフォルトではrange-based-forが適用できません. これは次に述べます．
ここは，しれっとクラスをクラスにコピーしています．デフォルトでは，こういう書き方をすると次の行動を行うことを利用しています．
- クラスをクラスにコピーするとは，メンバー変数の値をコピーすることである
- デフォルト挙動をあなたが不注意にプログラムしなおすと, 計算速度が落ちるだけなので, やめましょう．
他の場所は，まあ，計算しているだけですね．diffの最大値の求め方はパターンなので覚えましょう．

さて, この関数がちゃんと動作するには, myArray<double>クラスを強化して, range-based-forが動作するように改善しないとです．こんな感じです:

template <class T> class myArray
{
  ...
  std::vector<T> _data;
public:
  ...
  T& operator[](const int k)
  {
    return _data[k];
  }
  typename std::vector<T>::iterator begin()
  {
    return _data.begin();//内部配列 _data のbegin()を返却してお茶を濁す
  }
  typename std::vector<T>::iterator end()
  {
    return _data.end();//内部配列 _data のend()を返却してお茶を濁す
  }
  ...

このように, begin()とend()関数を持つクラスでは, range-based-forが動作します. 全部自作で作ると, いろいろと面倒なのですが, もともとstd::vector<T>に正しいbegin(), end()が定義されているので, そいつを転送するだけで十分です．

解答例2

JacobiのLaplaceは解く！