bj2_template

template.cpp

//C++11
#include <iostream>
#include <cmath>
#include <cstdio>
#include <vector>
#include <string>
#include <tuple>
#include <functional>

//備忘：問題によってacc函数の定義とmain函数内を書き換えれば良い、はず
using namespace std;

class V2 //2次元実数ベクトル
{
private:
	double x0,x1;
public:
	explicit V2 (double x0_=0.0, double x1_=0.0){
		x0=x0_;
		x1=x1_;
	}
	double& operator[](int i)//書き込み用添字演算子
	{
		switch(i){
		case 0:
			return x0;
		case 1:
			return x1;
		default:
			return x0;
		}
	}
	double operator[](int i) const//読み取り用添字演算子
	{
		switch(i){
		case 0:
			return x0;
		case 1:
			return x1;
		default:
			return x0;
		}
	}
	V2 operator+(const V2& rhm) const {//和
		return V2((*this)[0]+rhm[0],(*this)[1]+rhm[1]);
	}
	V2 operator-() const {//逆元
		return V2(-(*this)[0],-(*this)[1]);
	}
	double operator*(const V2& rhm) const {//内積
		return (*this)[0]*rhm[0]+(*this)[1]*rhm[1];
	}
	V2 operator*(double rhm) const {//スカラー倍
		return V2((*this)[0]*rhm,(*this)[1]*rhm);
	}
	double norm() const {//ノルム
		return std::sqrt((*this)*(*this));
	}
};

class V1 //1次元実数ベクトル
{
private:
	double x0;
public:
	explicit V1 (double x0_=0.0){
		x0=x0_;
	}
	double& operator[](int i)//書き込み用添字演算子。変ではあるけど多次元の場合にあわせて定義しておく
	{
		return x0;
	}
	double operator[](int i) const//読み取り用添字演算子
	{
		return x0;
	}
	V1 operator+(const V1& rhm) const {//和
		return V1((*this)[0]+rhm[0]);
	}
	V1 operator-() const {//逆元
		return V1(-(*this)[0]);
	}
	double operator*(const V1& rhm) const {//内積
		return (*this)[0]*rhm[0];
	}
	V1 operator*(double rhm) const {//スカラー倍
		return V1((*this)[0]*rhm);
	}
	double norm() const {//ノルム
		return abs((*this)[0]);
	}
};

class V3 //3次元実数ベクトル
{
private:
	double x0,x1,x2;
public:
	explicit V3 (double x0_=0.0, double x1_=0.0, double x2_=0.0){
		x0=x0_;
		x1=x1_;
		x2=x2_;
	}
	double& operator[](int i)//書き込み用添字演算子
	{
		switch(i){
		case 0:
			return x0;
		case 1:
			return x1;
		case 2:
			return x2;
		default:
			return x0;
		}
	}
	double operator[](int i) const//読み取り用添字演算子
	{
		switch(i){
		case 0:
			return x0;
		case 1:
			return x1;
		case 2:
			return x2;
		default:
			return x0;
		}
	}
	V3 operator+(const V3& rhm) const {//和
		return V3((*this)[0]+rhm[0],(*this)[1]+rhm[1],(*this)[2]+rhm[2]);
	}
	V3 operator-() const {//逆元
		return V3(-(*this)[0],-(*this)[1],-(*this)[2]);
	}
	double operator*(const V3& rhm) const {//内積
		return (*this)[0]*rhm[0]+(*this)[1]*rhm[1]+(*this)[2]*rhm[2];
	}
	V3 operator*(double rhm) const {//スカラー倍
		return V3((*this)[0]*rhm,(*this)[1]*rhm,(*this)[2]*rhm);
	}
	double norm() const {//ノルム
		return std::sqrt((*this)*(*this));
	}
	static V3 clpr(const V3& l, const V3& r){//外積
		const double x=l[1]*r[2]-l[2]*r[1];
		const double y=l[2]*r[0]-l[0]*r[2];
		const double z=l[0]*r[1]-l[1]*r[0];
		return V3(x,y,z);
	}
};

template <typename T>
class State //系の状態
{
public:
	T x; //座標
	T v; //速度
	State(T x_, T v_){
		x = x_;
		v = v_;
	}
};

const double NP = exp(1);

//加速度を求める
V3 acc(const State<V3> &s){
	const double r3 = pow(s.x.norm(),3);
	return V3::clpr(s.v,s.x)*(1/r3);
}

//オイラー法
template <typename T>
State<T> euler(const State<T> &s, double dt){
	V2 newx = s.x + s.v * dt;
	V2 newv = s.v + acc(s)*dt;
	return State<T>(newx,newv);
}

//2次ルンゲ・クッタ法(修正オイラー法)
template <typename T>
State<T> runge2(const State<T> &s, double dt){
	T newx1 = s.x + s.v * (dt/2);
	T newv1 = s.v + acc(s)*(dt/2);
	T newx2 = newx1 + newv1*(dt/2);
	T newv2 = newv1 + acc(State<T>(newx1,newv1))*(dt/2);
	return State<T>(newx2,newv2);
}

//4次ルンゲ・クッタ法(Rungeの原公式)
template <typename T>
State<T> runge4(const State<T> &s, double dt){
	T newx1 = s.x + s.v * (dt/2);
	T newv1 = s.v + acc(s)*(dt/2);
	T newx2 = newx1 + newv1*(dt/2);
	T newv2 = newv1 + acc(State<T>(newx1,newv1))*(dt/2);
	T newx3 = newx2 + newv2*dt;
	T newv3 = newv2 + acc(State<T>(newx2,newv2))*dt;
	T newx = (s.x+newx1*2.0+newx2*2.0+newx3)*(1/6.0);
	T newv = (s.v+newv1*2.0+newv2*2.0+newv3)*(1/6.0);
	return State<T>(newx,newv);
}

//1次シンプレクティック法（シンプレクティック・オイラー法）
template <typename T>
State<T> sym1(const State<T> &s, double dt){
	T newx = s.x + s.v * dt;
	T newv = s.v + acc(State<T>(newx,s.v))*dt;
	return State<T>(newx,newv);
}

//2次シンプレクティック法
template <typename T>
State<T> sym2(const State<T> &s, double dt){
	T newx = s.x + s.v * (dt/2);
	T newv = s.v + (-acc(State<T>(newx,s.v))*dt);
	T newx2 = newx + newv*(dt/2);
	return State<T>(newx2,newv);
}

//4次シンプレクティック法
template <typename T>
State<T> sym4(const State<T> &s, double dt){
	const double u=1.0/(2.0-cbrt(2.0));
	T newx1 = s.x + s.v * (u*dt/2.0);
	T newv1 = s.v + (-acc(State<T>(newx1,s.v))*u*dt);
	T newx2 = newx1 + newv1*(u*dt/2.0);
	T newx3 = newx2 + newv1*((1.0-2.0*u)*dt/2.0);
	T newv2 = newv1 + (-acc(State<T>(newx3,newv1))*(1.0-2.0*u)*dt);
	T newx4 = newx3 + newv2*((1.0-2.0*u)*dt/2.0);
	T newx5 = newx4 + newv2*(u*dt/2.0);
	T newv3 = newv2 + (-acc(State<T>(newx5,newv2))*u*dt);
	T newx6 = newx5 + newv3*(u*dt/2.0);
	return State<T>(newx6,newv3);
}


//時間発展（2次元）
void develop(State<V2> &state,double dt, double limit, function<State<V2>(const State<V2>&, double)> func){
	State<V2> s = state;
	double time=0;
	for( ; time <= limit; ){
		printf("%2.2f %e %e %e %e\n", time, s.x[0], s.x[1], s.v[0], s.v[1]);
		s = func(s,dt);
		time += dt;
	}
	printf("\n\n");
}

//時間発展（1次元）
void develop(State<V1> &state,double dt, double limit, function<State<V1>(const State<V1>&, double)> func){
	State<V1> s = state;
	double time=0;
	for( ; time <= limit; ){
		printf("%2.2f %e %e\n", time, s.x[0], s.v[0]);
		s = func(s,dt);
		time += dt;
	}
	printf("\n\n");
}
//時間発展（3次元）
void develop(State<V3> &state,double dt, double limit, function<State<V3>(const State<V3>&, double)> func){
	State<V3> s = state;
	double time=0;
	for( ; time <= limit; ){
		printf("%2.2f %e %e %e %e %e %e\n", time, s.x[0], s.x[1], s.x[2], s.v[0], s.v[1], s.v[2]);
		s = func(s,dt);
		time += dt;
	}
	printf("\n\n");
}

int main(void){
	double time=0; //時刻
	State<V3> state(V3(1.0,0.0,0.0), V3(-1.0,0.1,0.0)); //系の初期状態
	const double dt = 0.001; //刻み幅
	const double limit = 10.0;//範囲
	//develop(state,dt,limit,function<State<V3>(const State<V3>&, double)>(euler<V3>));
	//develop(state,dt,limit,function<State<V3>(const State<V3>&, double)>(runge2<V3>));
	//develop(state,dt,limit,function<State<V3>(const State<V3>&, double)>(runge4<V3>));
	//develop(state,dt,limit,function<State<V3>(const State<V3>&, double)>(sym1<V3>));
	//develop(state,dt,limit,function<State<V3>(const State<V3>&, double)>(sym2<V3>));
	develop(state,dt,limit,function<State<V3>(const State<V3>&, double)>(sym4<V3>));
	return 0;
}