#include <iostream>

#include <vector>

#include <cmath>

#include <fstream>

using namespace std;

vector<double> operator+ (vector<double>, vector<double>);

//vector<double> operator- (vector<double>, vector<double>);

vector<double> operator+= (vector<double> &, vector<double>);

vector<double> operator* (vector<double>, double);

ostream& operator« (ostream&, const vector<double>&);

vector<double> operator+ (vector<double> v, vector<double> v1)

{

return v += v1;

}

vector<double> operator- (vector<double> &v, vector<double> v1)

{

for(vector<double>::iterator iter = v.begin(), i = v1.begin(); iter != v.end(); iter++, i++)

{

(*iter) - (*i);

}

return v;

}

vector<double> operator+= (vector<double> &v, vector<double> v1)

{

for(vector<double>::iterator iter = v.begin(), i = v1.begin(); iter != v.end(); iter++, i++)

{

(*iter) += (*i);

}

return v;

}

vector<double> operator* (vector<double> v, double d)

{

for(vector<double>::iterator iter = v.begin(); iter < v.end(); iter++) *iter *= d;

return v;

}

ostream& operator« (ostream& os, const vector<double>& v)

{

os « "(";

for (vector<double>::const_iterator iter = v.begin(); iter != v.end() - 1; iter++)

{

os « *iter « "; ";

}

os « *(v.end() - 1) « ")";

return os;

}

vector<double> f(double t, vector<double> vect)

{

vector<double> vec(3);

vec[0] = vect[0]*2.0 + vect[1];

vec[1] = vect[0]+vect[1]*3.0-vect[2];

vec[2]=vect[1]*2.0+vect[2]*3.0-vect[0];

return vec;

}

void RungeCutt(vector<double> t, double h, int n, vector<vector<double> >& vect)

{

vector<vector<double> > K(4, vector<double>(3));

vector<vector<double> > temp(4, vector<double>(3));

for(int i = 0; i < n; i ++)

{

K[0] = f(t[i], vect[i])*h;

K[1] = f(t[i] + (h/2.0), vect[i] + K[0]*(1.0/2.0))*h;

K[2] = f(t[i] + (1.0/2.0)*h, vect[i] + K[1]*(1.0/2.0))*h;

K[3]= f(t[i]+h, vect[i]+K[2])*h;

vect[i + 1] = (K[0] + K[1]*2.0+K[2]*2.0+K[3])*(1.0/6.0) + vect[i];

}

}

vector<vector<double> > func(vector<double> t, const vector<double> vect)

{

int n = t.size();

double fs,sc,r,C1, C2,C3;

vector<vector<double> > vect1(n, vector<double>(3));

vect1[0] = vect;

fs = vect[0]+vect[2];

sc = vect[0]*3.0-vect[2]-vect[1];
r=vect[0]+vect[1];

C1 = sc*exp(t[0]*4.0)*(1.0/2.0);

C2 = (fs*(cos(t[0])+sin(t[0]))+vect[1]*(sin(t[0])-cos(t[0])))*exp(t[0]*3.0)*((-1.0)/2.0);

C3=(r*(sin(t[0])-cos(t[0]))-vect[2]*(cos(t[0])-sin(t[0])))*exp(t[0]*3.0)*((-1.0)/2.0);
cout«"C1="«C1«endl;
cout«"C2="«C2«endl;
cout«"C3="«C3«endl;

for (int i = 1 ; i < n; i++)

{

vect1[i][0] = C1*exp(t[i]*2.0) + C2*exp(t[i]*3.0)*cos(t[i])+ C3*exp(t[i]*3.0)*sin(t[i]);

vect1[i][1] = C2*exp(t[i]*3.0)*(cos(t[i])-sin(t[i]))+ C3*exp(t[i]*3.0)*(sin(t[i])+cos(t[i]));

vect1[i][2]=C1*exp(t[i]*2.0)+C2*exp(t[i]*3.0)*(cos(t[i])*2.0+sin(t[i]))+C3*exp(t[i]*3.0)*(sin(t[i])*2.0-cos(t[i]));

}

return vect1;

}

int main()

{

int n;

double a, b, h;

double Max = 0.0;

cout « "a = ";

cin » a;

cout « "b = ";

cin » b;

cout « "n = ";

cin » n;

h = (b - a)/n;

cout « "h = " « h « endl;

vector<vector<double> > vect(n + 1, vector<double>(3));

vector<vector<double> > vect1(n + 1, vector<double>(3));

vector<double> t(n + 1);

double k = a;

for (int i = 0; i <= n ; i++)

{

t[i] = k;

k += h;

}

cout « "x0 = ";

cin » vect[0][0];

cout « "y0 = ";

cin » vect[0][1];

cout«"z0=";

cin»vect[0][2];

RungeCutt(t, h, n, vect);

double p=0.0;

vect1 = func(t, vect[0]);

for(int i = 0; i <= n; i++)

{

cout «"t[" « i « "] = " « t[i] « endl;

cout « "u(x,y,z) = " « vect[i] « endl;

cout « "u1(x,y,z) = " « vect1[i] « endl;

if((p = fabs(vect[i][0] - exp(t[i]*2.0))) > Max ) Max = p;

if((p = fabs(vect[i][1] - 0 )) > Max) Max = p;

if((p=fabs( vect[i][2]-exp(t[i]*2.0)))>Max) Max=p;

cout « "u2(x,y,z) = " « exp(t[i]*2.0)«" "«0«" "«exp(t[i]*2.0) « endl;
}

cout « "\n(x(0), y(0),z(0)) = " « vect[0] « "\na = " « a « "\nb = " « b « "\nn = " « n « "\nmax eps = " « Max « endl;

return 0;

}