C++でルンゲクッタ法

ルンゲクッタ(RungeKutta)法のプログラム。

局所誤差を求めている。

#include<iostream>
#include<math.h>
 
using namespace std;
 
#define VLEN 3
#define REAL double
 
class vector{
private:
  REAL element[VLEN];
public:
  vector(REAL c = 0){
    for(int i = 0; i<VLEN; i++){
      element[i] = c;
    }
  }
  double & operator [] (int i){
    return element[i];
  }
 
  vector & operator += (const vector & b){
    for(int i = 0; i < VLEN; i++){
      element [i] += b.element[i];
    }
    return *this;
  }
 
  friend const vector operator + (const vector &, const vector &);
  friend const vector operator - (const vector &, const vector &);
  friend const vector operator * (const REAL &, const vector &);
  friend const double operator * (const vector &, const vector &);
  friend const vector operator & (const vector &, const vector &);
  friend ostream & operator << (ostream &, const vector &);
  friend istream & operator >> (istream &, vector &);
 
};
 
inline ostream & operator << (ostream & s, const vector & v){
  for(int i = 0; i < VLEN; i++){
    s << v.element[i] << " ";
  }
  return s;
}
 
inline istream & operator >> (istream & s, vector & v){
  for (int i = 0; i < VLEN; i++){
    s >> v.element[i];
  }
  return s;
}
 
inline const vector operator + (const vector &v1, const vector &v2){
  vector v3;
  for(int i = 0; i < VLEN; i++ ){
    v3.element[i] = v1.element[i] + v2.element[i];
  }
  return v3;
}
 
inline const vector operator - (const vector &v1, const vector &v2){
  vector v3;
  for(int i = 0; i < VLEN; i++ ){
    v3.element[i] = v1.element[i] - v2.element[i];
  }
  return v3;
}
 
inline const vector operator * (const REAL &c, const vector &v1){
  vector v3;
  for(int i = 0; i < VLEN; i++ ){
    v3.element[i] = v1.element[i] * c;
  }
  return v3;
}
 
inline const double operator * (const vector &v1, const vector &v2){
  double sum = 0;
  for(int i = 0; i < VLEN; i++ ){
    sum += v1.element[i] * v2.element[i];
  }
  return sum;
}
 
inline const vector operator & (const vector &v1, const vector &v2){
  vector v3;
  v3[0] = v1.element[1] * v2.element[2] - v1.element[0] * v2.element[1];
  v3[1] = v1.element[2] * v2.element[0] - v1.element[0] * v2.element[2];
  v3[2] = v1.element[0] * v2.element[1] - v1.element[1] * v2.element[0];
  return v3;
}
 
void rungekutta(vector, vector&, REAL, REAL);
void func(vector&, vector, REAL);
 
main(){
  int n;
  vector x, xnew, tmpx;
  REAL t, h, y;
  cout << "x = ";
  cin >> x;
  tmpx = x;
  for(float h = 0.000001; h <= 1; h *= 10){
    t = 0.0;
    x = tmpx;
    int j = 4/h;
    for(int i = 0; i < j; i++){
      rungekutta(x, xnew, t, h);
      x = xnew;
      t += h;
      if(i == 0){
        y = exp(t);
        cout << log10(h) << " " << t << " " << log10(sqrt((y - xnew[0])*(y - xnew[0]))) << endl;
      }
    }
  }
}
 
void rungekutta(vector x, vector &xnew, REAL t , REAL h){
  static vector u;
  static vector u2;
  static vector u3;
  static vector u4;
  func(u, x, t);
  func(u2, x+0.5*h*u, t+0.5*h);
  func(u3, x+0.5*h*u2, t+0.5*h);
  func(u4, x+h*u3, t+h);
  xnew = x + (1.0/6.0)*h*(u + 2*u2 + 2*u3 + u4);
}
 
void func(vector &u, vector x, REAL t){
  double a = 1;
  u = a * x;
}
 

ルンゲクッタのステップサイズの違いによる誤差の変化。