C++で中点法

授業で作った適当な中点法のプログラム。

使わないのにベクトルクラスを無駄に定義してる。
hがステップサイズで、xが初期値。

#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 euler(vector, vector&, REAL, REAL);
void func(vector&, vector, REAL);
 
main(){
  int n;
  vector x, xnew;
  REAL t, h;
  cout << "x = ";
  cin >> x;
  cout << "h = ";
  cin >> h;
  t = 0.0;
  for(int i = 0; i < 4/h; i++){
    euler(x, xnew, t, h);
    cout << t << " " << x << endl;
    x = xnew;
    t += h;
  }
  cout << t << " " << xnew << endl;
}
 
void euler(vector x, vector &xnew, REAL t , REAL h){
  static vector u;
  static vector u2;
  func(u, x, t);
  func(u2, x+h*u, t+h);
  xnew = x + 0.5*h*u + 0.5*h*u2;
}
 
void func(vector &u, vector x, REAL t){
  double a = 1;
  u = a * x;
}
 

ステップサイズによる変化は下のような感じになった。