Добавил:

Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.

Вуз:

Уральский Федеральный университет им. Б.Н. Ельцина «УПИ»

Предмет:

[НЕСОРТИРОВАННОЕ]

Файл:

Diplom_SamoylovAL.pdf

Скачиваний:

Добавлен:

23.02.2015

Размер:

7.89 Mб

Скачать

☆

<<< < Предыдущая 1 2 3 4 5 6 7 89 / 99

6 ЗАКЛЮЧЕНИЕ

Разработанная программа может применяться для решения достаточно ши рокого класса задач управления, а также задачи Коши для уравнения Гамиль тона-Якоби-Беллмана для случаев фазового пространства 1, 2, 3. Предпола гается выполненными условия существования, единственности и продолжимости классических характеристик.

Изучена библиотека параллельных вычислений OpenMP. Распараллелены ал горитмы попятного интегрирования и поиска пересечений характеристик.

Получены численные решения ряда модельных задач. Для задач построены оптимальные процессы, функция цены и сеточный оптимальный синтез.

Возможными направлениями развития работы являются повышение размер ности решаемых задач, вопросы оптимального подбора параметров, а также оп тимизация программы по времени счета.

Поставленная задача решена полностью.

Благодарности

Работа выполнена при финансовой поддержке Российского фонда фундамен тальных исследований (проект №11–01–00214).

7 СПИСОК ИСПОЛЬЗОВАННЫХ

ИСТОЧНИКОВ

1.Subbotina N. N. The method of characteristics for Hamilton–Jacobi equations and applications to dynamical optimization. NY: Springer, 2006. Vol. 135. Pp. 2955–3091.

2.www.gnu.org/software/gsl.

3.www.cgal.org.

4.www.openmp.org.

5.Понтрягин Л. С., Болтянский В. Г., Гамкрелидзе Р. В., Мищенко Е. Ф. Мате матическая теория оптимальных процессов. Москва: Наука, 1961.

6.Bellman R. Dynamic Programming. Princeton: Univ. Press, 1957.

7.Красовский Н. Н. Теория управления движением. Москва: Наука, 1968.

8.Красовский Н. Н., Субботин А. И. Позиционные дифференциальные игры. Москва: Наука, 1974.

9.Петровский И. Г. Лекции по теории обыкновенных дифференциальных урав нений. Москва: Изд–во МГУ, 1984.

10.Курант Р. Уравнения с частными производными. Москва: Мир, 1964.

11.Clarke F. H., Vinter R. The relationship between the maximum principle and dynamic programming // S.I.A.M. J. Contr. Optimiz. 1987. Vol. 5. Pp. 1291–1311.

12.Barron E. N., Jensen R. The Pontryagin maximum principle from dynamic pro gramming and viscosity solutions to first–order partial differential equations // Trans. Amer. Math. Soc. 1986. Vol. 298, no. 2. Pp. 635–641.

13.Mirica S. Extending Cauchy’s method of characteristics for Hamilton–Jacobi equa tions // Stud. Cerc. Mat. 1985. Vol. 37, no. 6. Pp. 555–565.

14.Субботин А. И. Обобщенные решения уравнений в частных производных пер вого порядка. Перспективы динамической оптимизации. Москва; Ижевск: Ин ститут компьютерных исследований, 2003.

15.Crandall M. G., Lions P. L. Viscosity solutions of Hamiltion–Jacobi equations // Trans. Amer. Math. Soc. 1983. Vol. 277. Pp. 1–42.

16.Препарата Ф., Шеймос М. Вычислительная геометрия: Введение. Москва: Мир, 1989.

17.Crandall M. G., Lions P.-L. Viscosity Solutions of Hamilton–Jacobi Equations // Trans. Am. Math. Soc. 1983. Vol. 277:1. Pp. 1–42.

18.Филиппов А. Ф. О некоторых вопросах теории оптимального регулирова ния // Вестник Моск. Ун.-та. 1959. № 2. С. 25–32.

19.Болтянский В. Г. Математические методы оптимального управления. Москва: Наука, 1969.

20.Пацко В. С. Поверхности переключения в линейных дифференциальных иг рах. Екатеринбург: ИММ УрО РАН, 2004.

21.Субботина Н. Н. Обобщенный метод характеристик в задаче оптимального управления с липшицевыми входными данными // Изв. УрГУ: Математика и механика. 2006. Т. 4. С. 171–176.

22.Токманцев Т. Б. Обобщенный метод характеристик в решении задач опти мального управления с фиксированным моментом окончания: Кандидатская диссертация / ИММ УрО РАН. Екатеринбург, 2011.

23.Шолохович Ф. А. Лекции по дифференциальным уравнениям. Екатеринбург: Уральское изд–во, 2005.

24.Варга Д. Оптимальное управление дифференциальными и функциональными уравнениями. Москва: Наука, 1977.

25.Филиппов А. Ф. Дифференциальные уравнения с разрывной правой частью. Москва: Наука, 1985.

8 ПРИЛОЖЕНИЕ

//MAIN.CPP

#include "Problem2.h" #include "TSynthes2.h" #include <iostream> #include <time.h>

int main()

{

std::cout <<"Problem - "<<ProblemName<<std::endl <<"Np = "<<NP<<std::endl;

time_t start, end; time(&start);

TBuilder2 :: TPointer VF = Construct(); std::ostringstream DirName;

DirName << "var_" <<ProblemName; system((string("mkdir ") + DirName.str()).c_str()); VF->Instant2Files(DirName.str());

time(&end);

double Time_Clock = difftime(end,start); std::cout<<Time_Clock<<std::endl

<<"Construct completed"<<std::endl; time(&start);

TSynthes2 S(VF); S.BuildTrajectories(SynthesGrid); S.Instant2Files(DirName.str()); S.Instant2MatlabFile(DirName.str()); time(&end);

Time_Clock = difftime(end,start); std::cout<<Time_Clock<<std::endl

<<"Synthes completed"<<std::endl;

return 0;

}

//PROBLEM2.H #pragma once #include <math.h> #include "TBuilder2.h" #include "TFactory2.h"

#include "Definitions2.h"

using namespace std;

namespace SpherPendulum{

extern std::string	ProblemName;
const double	Eps = 1e-4;
const double	dt = 0.05;
const double	dx = 0.2;//dt*dt;
const double	rho = 0;
const double	T = 1;
//extern double	U;
const Vector2	x1(0,-5),
	x2(0,5);

vector<Vector2> CreateGrid(); extern vector<Vector2> SynthesGrid;

const int XNum = static_cast<int>((x2.x() - x1.x()) / dx); const int TNum = static_cast<int>(T / dt);

double Sigma(const Vector2 &); Vector2 DxSigma(const Vector2 &); int RHS(

double tau, const double y[], double result[], void *params

);

double Control(const double t, const Vector2 &x, const Vector2 &p);

Vector2 SynthesControl(const double t, const Vector2 &Point, const Vector2 &Target);

Vector2 f(const double t, const Vector2 &x, const Vector2 &u);//вторая координата u фиктивна double g(const double t, const Vector2 &x, const Vector2 &u);//вторая координата u фиктивна

}
namespace Didona{
extern std::string	ProblemName;
const double	Eps = 1e-4;
const double	dt = 0.1;
const double	dx = 0.5;//dt*dt;
const double	rho = 0;
const double	T = 1;
//extern double	U;
const Vector2	x1(-1,-1),
	x2(1,1);

vector<Vector2> CreateGrid(); extern vector<Vector2> SynthesGrid;

const int XNum = static_cast<int>((x2.x() - x1.x()) / dx); const int TNum = static_cast<int>((T+Eps) / dt);

double Sigma(const Vector2 &); Vector2 DxSigma(const Vector2 &); int RHS(

double tau, const double y[], double result[], void *params

);

Vector2 Control(const double t, const Vector2 &x, const Vector2 &p);

Vector2 SynthesControl(const double t, const Vector2 &Point, const Vector2 &Target); Vector2 f(const double t, const Vector2 &x, const Vector2 &u);

double g(const double t, const Vector2 &x, const Vector2 &u);

}

namespace GravityProblem{

static double v_T = -98.0;//!!!

extern std::string	ProblemName;
const double	Eps = 1e-4;
const double	dt = 0.1;
const double	dx = dt*dt;
const double	rho = 0;
const double	T = 10;
extern double	U;
const Vector2	x1(-1e-1, v_T - 1e-1),

x2( 1e-1, v_T + 1e-1); const int XNum = static_cast<int>((x2.x() - x1.x()) / dx); const int TNum = static_cast<int>(T / dt);

double y(const double t);

double Integrant(const double t, const Vector2 &x, const double u); double Control(const double t, const Vector2 &x, const Vector2 &p);

Vector2 f(const double t, const Vector2 &x, const Vector2 &p, const double u); double g(const double t, const Vector2 &x, const Vector2 &p, const double u); const Vector2 h(const double t, const Vector2 &x, const Vector2 &p);

double Sigma(const Vector2 &); Vector2 DxSigma(const Vector2 &); int RHS(

double tau, const double y[], double result[], void *params

);

}

using namespace SpherPendulum;

TBuilder2 :: TPointer Construct();

//PROBLEM2.CPP

#include"Problem2.h"

#define _USE_MATH_DEFINES #include "math.h"

double sqr(const double x){ return x*x;

}

namespace SpherPendulum

{

vector<Vector2> SynthesGrid; vector<Vector2> CreateGrid()

{

vector<Vector2> result; for(int j = 0; j <= XNum; ++j)

{

const Vector2 y(0, x1.y() + j * dx); result.push_back(y);

}

return result;

}

std::string ProblemName = "Pendulum"; const double grav = 9.8;

double Control(const double t, const Vector2 &x, const Vector2 &p)

{

double u0 = -p.y()*sin(x.x())/Eps; if(abs(u0) <= 1)

return u0; return u0/abs(u0);

}

Vector2 SynthesControl(const double t, const Vector2 &Point, const Vector2 &Target)

{

Vector2 p = Target - Point; if(p.y()*sin(Point.x()) < 0)

return Vector2(-1,0); //вторая координата фиктивна

else

return Vector2(1,0);

}

Vector2 f(const double t, const Vector2 &x, const Vector2 &u)//вторая координата u фиктивна

{

return Vector2(x.y(),(-grav+u.x())*sin(x.x())+2*sin(2*x.x()));

}

double g(const double t, const Vector2 &x, const Vector2 &u)//вторая координата u фиктивна

{

return Eps*u.x()*u.x()/2;

}

const Vector2 h(const double t, const Vector2 &x, const Vector2 &p, const double &u)

{

return Vector2(p.y()*((-grav+u)*cos(x.x())+4*cos(2*x.x())), p.x());

}

double Sigma(const Vector2 &x)

{

return x.x()*x.x()/2;

}

Vector2 DxSigma(const Vector2 &x)

{

return Vector2(x.x(),0);

}

int RHS( double tau,

const double y[], double result[], void *params

){

double t = T - tau;

Vector2 x = Vector2(y[0], y[1]); Vector2 p = Vector2(y[2], y[3]); const double u = Control(t,x,p); Vector2 X = f(t,x,Vector2(u,0)); Vector2 P = h(t,x,p,u);

result[0] = -X.x(); result[1] = -X.y(); result[2] = P.x(); result[3] = P.y();

result[4] = g(t,x,Vector2(u,0)); return GSL_SUCCESS;

}

namespace Didona

{

vector<Vector2> SynthesGrid; vector<Vector2> CreateGrid()

{

vector<Vector2> result; for(int i = 0; i <= XNum; ++i)

{

for(int j = 0; j <= XNum; ++j)

{

const Vector2 y(x1.x() + i * dx, x1.y() + j * dx); result.push_back(y);

}

return result;

}

std::string ProblemName = "Didona";

Vector2 Control(const double t, const Vector2 &x, const Vector2 &p)

{

Vector2 u0(-(p.x()-x.y())/Eps,-(p.y()+x.x())/Eps); double l2 = u0.squared_length();

if(l2 <= 1)

return u0; return u0/sqrt(l2);

}

Vector2 SynthesControl(const double t, const Vector2 &Point, const Vector2 &Target)

{

double l = sqrt((Target-Point).squared_length()); return (Target-Point)/l;

}

Vector2 f(const double t, const Vector2 &x, const Vector2 &u)

{

return u;

}

double g(const double t, const Vector2 &x, const Vector2 &u)

{

return (x.x()*u.y() - x.y()*u.x()) + Eps*u.squared_length()/2;

}

const Vector2 h(const double t, const Vector2 &x, const Vector2 &p, const Vector2 &u)

{

return Vector2(u.y(), -u.x());

}

double Sigma(const Vector2 &)

{

return 0;

}

Vector2 DxSigma(const Vector2 &x)

{

return Vector2(0,0);

}

int RHS( double tau,

const double y[], double result[], void *params

){

double t = T - tau;

Vector2 x = Vector2(y[0], y[1]); Vector2 p = Vector2(y[2], y[3]); const Vector2 u = Control(t,x,p); Vector2 X = f(t,x,u);

Vector2 P = h(t,x,p,u); result[0] = -X.x(); result[1] = -X.y(); result[2] = P.x(); result[3] = P.y(); result[4] = g(t,x,u); return GSL_SUCCESS;

}

namespace GravityProblem

{

std::string ProblemName = "GravityProblem";

double U = 10;

double y(const double t)

{

const double g = 9.8;

const double h0 = g*10*10/2; return -g*t*t/2 + h0;

}

double Control(const double t, const Vector2 &x, const Vector2 &p)

{

double u0 = -p.y()/Eps; if(fabs(u0) < U)

return u0;

double H1 = p.x()*x.y() + p.y()*U + Integrant(t,x,U); double H2 = p.x()*x.y() + p.y()*(-U) + Integrant(t,x,-U); return H1 < H2 ? U:-U;

}

Vector2 f(const double t, const Vector2 &x, const Vector2 &p, const double u)

{

return Vector2(x.y(),u);

}

double Integrant(const double t, const Vector2 &x, const double u)

{

return sqr(x.x() - y(t))/2 + u*u*Eps/2;

}

double g(const double t, const Vector2 &x, const Vector2 &p, const double u)

{

return Integrant(t,x,u);

}

const Vector2 h(const double t, const Vector2 &x, const Vector2 &p)

{

return Vector2(x.x() - y(t), p.x());

}

double Sigma(const Vector2 &)

{

return 0;

}

Vector2 DxSigma(const Vector2 &x)

{

return Vector2(0,0);

}

int RHS(

double tau, const double y[], double result[], void *params

)

{

double t = T - tau;

Vector2 x = Vector2(y[0], y[1]); Vector2 p = Vector2(y[2], y[3]); const double u = Control(t,x,p); Vector2 X = f(t,x,p,u);

Vector2 P = h(t,x,p); result[0] = -X.x(); result[1] = -X.y(); result[2] = P.x(); result[3] = P.y(); result[4] = g(t,x,p,u); return GSL_SUCCESS;

}

TBuilder2 :: TPointer Construct()

{

vector<Vector2> Grid = CreateGrid(); SynthesGrid.clear();

for(int i = 0; i <= 2*XNum; ++i)

{

for(int j = 0; j <= 2*XNum; ++j)

{

const Vector2 y(x1.x() + i * dx/2, x1.y() + j * dx/2); SynthesGrid.push_back(y);

}

TFactory2 Factory(dt);

TSection2 Section = Factory.Create(T, Grid, Sigma, DxSigma, rho); TBuilder2 Builder(Section);

std::cout<<"0 ";

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

{

std::cout<<i<<endl<<" ";

Section = Factory.Create(Section); Builder.AddSection(Section);

}

return Builder.Get();

}

// TSynthes2.h #pragma once #include <assert.h>

#include "TValueFunction2.h" #include "Problem2.h" #include "TBuilder2.h"

#include <CGAL/Delaunay_triangulation_2.h>

#include <CGAL/Triangulation_vertex_base_with_info_2.h>

typedef CGAL::Triangulation_vertex_base_with_info_2<size_t, K> Vb;

typedef CGAL::Triangulation_data_structure_2<Vb>	Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds>	Triangulation;
typedef Triangulation::Point Point;

class TSynthes2

{

public:

TSynthes2(const TBuilder2 :: TPointer VF); void BuildTriang();

Vector2 GetNearestPoint(size_t t, const Vector2 &Point);

Vector2 Aim(size_t t, const Vector2 &Point, const Vector2 &Target, Vector2 &u); void BuildTrajectory(Vector2 Point);

void BuildTrajectories(const vector<Vector2> &Array); void Instant2Files(const string &DirName);

void Instant2MatlabFile(const string &DirName);

private:

vector<Triangulation> Triangs; TBuilder2 :: TPointer VF; vector<vector<Vector2> > Trajectories; vector<vector<Vector2> > UControls; vector<double> Values;

};

//TSYNTHES2.CPP #include "TSynthes2.h" #include <omp.h>

void TSynthes2::BuildTriang()

{

#pragma omp parallel for num_threads(NP) for(int t = 0; t < TNum; ++t)

{

vector<Vector2> Array = VF->GetSectionPoints(t); std::vector< std::pair<Point,size_t> > points; for(size_t i=0; i < Array.size();++i)

points.push_back( std::make_pair(CGAL::ORIGIN+Array[i],i) ); Triangulation Triang;

Triang.insert( points.begin(),points.end() ); Triangs[t] = Triang;

}

TSynthes2::TSynthes2(const TBuilder2 :: TPointer _VF) :VF(_VF),Triangs(TNum){}

Vector2 TSynthes2::GetNearestPoint(size_t t, const Vector2 &Point)

{

Triangulation::Vertex_handle Vertex = Triangs[t].nearest_vertex(CGAL::ORIGIN+Point); return VF->GetSectionPoints(t+1)[Vertex->info()];

}

Vector2 TSynthes2::Aim(size_t t, const Vector2 &Point, const Vector2 &Target, Vector2 &u)

{

assert(t < TNum); if(Point == Target)

u = Vector2(0,0);

else

u = SynthesControl(t*dt, Point, Target); return f(t*dt, Point, u);

}

void TSynthes2::BuildTrajectory(Vector2 Point)

{

Vector2 Target; Vector2 Velocity;

vector<Vector2> Trajectory; vector<Vector2> UControl; double Value=0;

Vector2 u; Trajectory.push_back(Point); for(size_t i = 0; i < TNum; ++i)

{

Target = GetNearestPoint(i,Point); Velocity = Aim(i, Point, Target, u); UControl.push_back(u); Value+=g(i*dt,Point,u);

Point = Point + Velocity * dt; Trajectory.push_back(Point);

}

Value-=g(0,Trajectory[0],UControl[0])/2; Value-=g(T,Trajectory[TNum],UControl[TNum-1])/2; Value*=dt;

Value+=Sigma(Point); Values.push_back(Value); Trajectories.push_back(Trajectory); UControl.push_back(Vector2(0,0)); UControls.push_back(UControl);

}

void TSynthes2::BuildTrajectories(const vector<Vector2> &Array)

{

BuildTriang();

for(size_t i = 0; i < Array.size(); ++i)

{

std::cout<<i<<endl;

BuildTrajectory(Array[i]);

}

void TSynthes2::Instant2Files(const string &DirName)

{

assert(!Trajectories.empty());

unsigned int Size = Trajectories[0].size(); for(unsigned int t = 0; t < Size; ++t){

ostringstream FileName;

FileName << DirName+"/s_" << t << ".txt"; ofstream Out(FileName.str().c_str());

Out << 3 << " " << Trajectories.size() << endl;

for(unsigned int k = 0; k < Trajectories.size(); ++k){ Out << Trajectories[k][t].x() << " ";

Out << Trajectories[k][t].y() << " "; Out << UControls[k][t].x() << " "; Out << UControls[k][t].y() << " "; Out << Values[k] << endl;

}

Out.close();

}

void TSynthes2::Instant2MatlabFile(const string &DirName)

{

ostringstream FileName;

FileName << DirName<<"_Synthes.m"; ofstream MatlabFile(FileName.str().c_str());

unsigned int Size = Trajectories[0].size();
MatlabFile << "X1 = zeros(" << Trajectories.size()		<< ", " << Size << ");\n";
MatlabFile << "X2 = zeros(" << Trajectories.size()		<< ", " << Size << ");\n";
MatlabFile << "U1 = zeros(" << UControls.size()		<< ", " << Size << ");\n";
MatlabFile << "U2 = zeros(" << UControls.size()		<< ", " << Size << ");\n";
MatlabFile << "Z	= zeros(" << Trajectories.size()	<< ", 1);\n";
MatlabFile << "for t=1:" << Size << ";\n";
MatlabFile << " Name = strcat('"+DirName+"/s_',num2str(t-1));\n";
MatlabFile << " Name = strcat(Name,'.txt');\n";
MatlabFile << "	File = fopen(Name, 'r');\n";
MatlabFile << "	fscanf(File,'%d',1);\n";
MatlabFile << " N = fscanf(File,'%d',1);\n";
MatlabFile << "	for i=1:N;\n";
MatlabFile << "	X1(i,t) = fscanf(File,'%g',1);\n";
MatlabFile << "	X2(i,t) = fscanf(File,'%g',1);\n";
MatlabFile << "	U1(i,t) = fscanf(File,'%g',1);\n";
MatlabFile << "	U2(i,t) = fscanf(File,'%g',1);\n";
MatlabFile << "	Z(i) = fscanf(File,'%g',1);\n";
MatlabFile << "	end\n";
MatlabFile << "	fclose(File);\n";

MatlabFile << "end\n";

MatlabFile << "figure(1);\n";

MatlabFile << "Time = 0 : "<<dt<<" : "<<T<<";\n";

MatlabFile << "hold on;\n";

MatlabFile << "for i=1:N;\n";

MatlabFile << "\tplot3(Time, X1(i,:), X2(i,:));\n";

MatlabFile << "end\n";

MatlabFile << "xlabel('t');\n";

MatlabFile << "ylabel('x_1');\n";

MatlabFile << "zlabel('x_2');\n";

MatlabFile << "figure(2);\n";

MatlabFile << "Time = 0 : "<<dt<<" : "<<T-dt<<";\n"; MatlabFile << "hold on;\n";

MatlabFile << "for i=1:N;\n";

MatlabFile << "\tplot3(Time, U1(i,1:"<<Size-1<<"), U2(i,1:"<<Size-1<<"));\n"; MatlabFile << "end\n";

MatlabFile << "xlabel('t');\n"; MatlabFile << "ylabel('u_1');\n"; MatlabFile << "zlabel('u_2');\n";

MatlabFile << "figure(3);\n";

MatlabFile << "Triang = delaunay(X1(:,1),X2(:,1));\n";

MatlabFile << "trimesh(Triang,X1(:,1),X2(:,1),Z(:));\n"; MatlabFile << "xlabel('x_1');\n";

MatlabFile << "ylabel('x_2');\n"; MatlabFile << "zlabel('S');\n";

}

// TBuilder2.h #pragma once

#include "TValueFunction2.h" #include "Tree2.h"

#include "../Include/Tree.h" #include "TFactory2.h"

#include <boost/tr1/memory.hpp>

class TBuilder2

{

public:

typedef std::tr1::shared_ptr<TValueFunction2> TPointer; TBuilder2(TSection2 &Section);

TPointer Get();

void AddSection(TSection2 &);

private:

bool FirstTime; TTree Times;

vector<TTree2> Trajectories; vector<TTree> Values; TPointer Ptr;

};

//TBUILDER2.CPP #include "TBuilder2.h"

TBuilder2::TBuilder2(TSection2 &Section) : Times(Section.Time())

{

const vector<Vector2> &X = Section.GetX(); const vector<double> &Z = Section.GetZ(); for(unsigned int i = 0; i < X.size(); ++i)

{

Trajectories.push_back(TTree2(X[i])); Values.push_back(TTree(Z[i]));

}

void TBuilder2::AddSection(TSection2 &Section)

{

Section.Formula();

const vector<Vector2> &X = Section.GetX(); const vector<double> &Z = Section.GetZ(); Times.Grow(Section.Time());

for(unsigned int i = 0; i < X.size(); ++i){ Trajectories[i].Grow(X[i]); Values[i].Grow(Z[i]);

}

TBuilder2::TPointer TBuilder2::Get()

{

TValueFunction2 * tmp = new TValueFunction2( Times,

Trajectories, Values

);

Ptr = TPointer(tmp); return Ptr;

}

//TFACTORY2.H #pragma once

#include "Definitions2.h" #include "TSection2.h" #include "TValueFunction2.h" #include <gsl/gsl_odeiv.h> #include <gsl/gsl_errno.h> #include <gsl/gsl_matrix.h> #include <vector>

using std :: vector;

class TFactory2

{

public:

TFactory2(double _dt); virtual ~TFactory2(void);

TSection2 Create(const TSection2&); TSection2 Create(

const double t,

const vector<Vector2> &X, TBoundaryFunction2 Sigma, TGradientFunction2 DxSigma, const double rho

);

private:

double dt;

};

//TFACTORY2.CPP #include "TFactory2.h" #include "Problem2.h" #include <assert.h> #include <algorithm> using std::transform;

TFactory2::TFactory2(double _dt) :dt(_dt)

{

}

TFactory2::~TFactory2(void){

}

void Integrate(

const double tau1, const double tau2, vector<double> &Y

){

assert(tau1 < tau2);

int Dimension = Y.size(); double *y = &Y[0]; double y_err[16],

dydt_in[16], dydt_out[16];

gsl_odeiv_system sys = { RHS, NULL, Dimension};

const gsl_odeiv_step_type *T = gsl_odeiv_step_rk8pd; gsl_odeiv_step * const s = gsl_odeiv_step_alloc(T, Dimension);

GSL_ODEIV_FN_EVAL (&sys, tau1, y, dydt_in);

const double h = 2E-5; double tau = tau1; while (tau < tau2){

int status = gsl_odeiv_step_apply (

s, tau, h, y, y_err, dydt_in, dydt_out, &sys

);

assert (status == GSL_SUCCESS); for(size_t i = 0; i < Dimension; ++i)

{

dydt_in[i] = dydt_out[i];

}

tau += h;

}

gsl_odeiv_step_free(s);

}

#include <time.h> #include <omp.h>

TSection2 TFactory2::Create(const TSection2 &Section){ vector<Vector2> X = Section.GetX(),

P = Section.GetP(); vector<double> Z = Section.GetZ();

double t2 = Section.Time(), t1 = t2 - dt;

int Size = static_cast<int>(X.size()); #pragma omp parallel for num_threads(NP) for(int i = 0; i < Size; ++i)// omp требует int

{

assert(omp_get_num_threads() <= omp_get_num_procs()); assert(omp_get_num_threads() == NP);

vector<double> Y(5); Y[0] = X[i].x();

Y[1] = X[i].y();

Y[2] = P[i].x();

Y[3] = P[i].y(); Y[4] = Z[i];

Integrate(T-t2, T-t1, Y); X[i] = Vector2(Y[0], Y[1]); P[i] = Vector2(Y[2], Y[3]); Z[i] = Y[4];

}

return TSection2(t1, X, P, Z, Section.GetRho());

}

TSection2 TFactory2::Create( const double t,

const vector<Vector2> &X, TBoundaryFunction2 Sigma, TGradientFunction2 DxSigma, const double rho

){

vector<Vector2> P; vector<double> Z;

transform(X.begin(), X.end(), back_inserter(P), DxSigma); transform(X.begin(), X.end(), back_inserter(Z), Sigma); return TSection2(t, X, P, Z, rho);

}

100

//DEFINITIONS.H #pragma once

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h> #include <vector>

using std :: vector; using std :: pair;

struct K : CGAL::Exact_predicates_inexact_constructions_kernel{}; typedef K::Vector_2 Vector2;

typedef Vector2(*TGradientFunction2)(const Vector2&); typedef double (*TBoundaryFunction2)(const Vector2&);

//typedef Vector2(*TFunction2)(const Vector2&, const Vector2&); typedef pair< vector<size_t>::iterator, vector<size_t>::iterator > TRange;

//TREE.H #pragma once #include <deque>

#include <iostream> #include <iomanip> #include <algorithm> #include <string> using std :: deque; using std :: istream; using std :: ostream; using std :: ofstream; using std :: endl; using std :: string;

class TTree{ public:

TTree(const double Root)

:_Branches(1, Root), _Root(_Branches.back()), _Stop(false){

}

TTree(const TTree &Tree) :_Root(Tree._Root),_Branches(Tree._Branches), _Stop(false){

}

TTree(istream &In);

TTree operator = (const TTree &Tree); void Save(ostream &Out) const;

void SaveTo(const string &FileName)const; bool operator <(const TTree &Tree) const; bool operator ==(const TTree &Tree) const; bool operator !=(const TTree &Tree) const;

const double& operator [](const size_t i) const{return _Branches[i];} bool IsCorrect();

bool Empty(){return _Branches.empty();} size_t Size(){return _Branches.size();} void Grow(const double);

void Stop();

private:

deque<double> _Branches; double _Root;

bool _Stop;

};

ostream& operator << (ostream &Out, const TTree &Tree);

//TREE.CPP

#include"Tree.h"

101

#include<assert.h>

#include<string>

#include<iostream>

#include<fstream> using std :: ofstream; using std :: string;

TTree TTree::operator = (const TTree &Tree){

_Branches = Tree._Branches; _Root = _Branches.back(); _Stop = Tree._Stop;

return *this;

}

void TTree::Save(ostream &Out) const{ //Out << "Trajectory = [" << endl;

for(unsigned int i = 0; i < _Branches.size(); ++i){

Out << ' ' << std :: fixed << std :: setprecision(20) << _Branches[i];

}

bool TTree::operator <(const TTree &Tree) const{ return _Root < Tree._Root;

}

void TTree::Grow(const double Branch){ if(!_Stop){

_Branches.push_front(Branch);

}

void TTree::Stop(){ _Stop = true;

}

bool TTree::operator ==(const TTree &Tree) const{

return !(_Root < Tree._Root) && !(Tree._Root < _Root) && (_Branches.size() == Tree._Branches.size());

}

bool TTree::operator !=(const TTree &Tree) const{ return !(*this == Tree);

}

bool TTree::IsCorrect(){

//?=> Уникальность корней return _Root == _Branches.back();

}

const deque<double> Split(const string &S){ std :: deque<double> Result;

std :: string::size_type PrevPos = 0;

while(true){

std :: string::size_type Pos = S.find(' ', PrevPos);

std :: string Number = S.substr(PrevPos, Pos - PrevPos);

PrevPos = Pos + 1; Result.push_back(atof(Number.c_str())); if(Pos == std :: string :: npos){

break;

}

return Result;

}

ostream& operator << (ostream &Out, const TTree &Tree){ Tree.Save(Out);

return Out;

102

}

TTree::TTree(istream &In)

{

std :: string Branches; In.ignore(100, '['); In.ignore(); //' ' getline(In, Branches);

assert(!Branches.empty()); _Branches = Split(Branches); _Root = _Branches.back();

}

void TTree::SaveTo(const string &FileName)const{

ofstream File(FileName.c_str(), std::ios::out|std::ios::binary ); double *Buf = new double [_Branches.size()]; copy(_Branches.begin(), _Branches.end(), Buf);

const char * Buffer = reinterpret_cast<const char*>(&Buf[0]); unsigned int Size = _Branches.size() * sizeof(Buf[0]); File.write(Buffer, Size);

File.close(); delete [] Buf;

}

//TREE2.H #pragma once #include <deque> #include <iostream> #include <iomanip>

#include <algorithm> #include <string> #include <vector> #include "Definitions2.h" using std :: deque;

using std :: istream; using std :: ostream; using std :: ofstream; using std :: endl; using std :: string;

class TTree2

{

public:

TTree2(const Vector2& Root) :_Branches(1, Root){

}

TTree2(const TTree2 &Tree) :_Branches(Tree._Branches){

}

TTree2 operator = (const TTree2 &Tree); void Save(ostream &Out) const;

void SaveTo(const string &FileName)const; bool operator ==(const TTree2 &Tree) const; bool operator !=(const TTree2 &Tree) const;

const Vector2& operator [](const size_t i) const{return _Branches[i];} void Grow(const Vector2&);

inline unsigned int Size(){return _Branches.size();}

friend ostream& operator << (ostream &Out, const TTree2 &Tree); Vector2 GetBranche(size_t i) const {return _Branches[i];}

private:

deque<Vector2> _Branches;

};

//TREE2.CPP

103

#include "Tree2.h" #include <assert.h> #include <string> #include <iostream> #include <fstream> using std :: ofstream; using std :: string;

using std :: ostream_iterator;

TTree2 TTree2::operator = (const TTree2 &Tree){ _Branches = Tree._Branches;

return *this;

}

void TTree2::Save(ostream &Out) const{ //Out << "Trajectory = [" << endl;

for(unsigned int i = 0; i < _Branches.size(); ++i){

Out << ' ' << std :: fixed << std :: setprecision(20) << _Branches[i];

}

void TTree2::Grow(const Vector2& Branch){ _Branches.push_front(Branch);

}

bool TTree2::operator ==(const TTree2 &Tree) const

{

return _Branches.size() == Tree._Branches.size();

}

bool TTree2::operator !=(const TTree2 &Tree) const

{

return !(*this == Tree);

}

ostream& operator << (ostream &Out, const TTree2 &Tree){ Tree.Save(Out);

return Out;

}

void TTree2::SaveTo(const string &FileName)const{ const size_t Size = 2*_Branches.size();

ofstream File(FileName.c_str(), std::ios::out|std::ios::binary ); double *Buf = new double [Size];

size_t i = 0;

for(deque<Vector2>::const_iterator it = _Branches.begin(); it != _Branches.end(); ++it)

{

Buf[i++] = it->x(); Buf[i++] = it->y();

}

const char * Buffer = reinterpret_cast<const char*>(&Buf[0]); File.write(Buffer, Size * sizeof(Buf[0]));

File.close(); delete [] Buf;

}

//TSECTION2.H #pragma once #include<vector> #include<iostream> #include "Definitions2.h"

#include "../Include/Functions.h" #include <time.h>

#include <omp.h> using std :: vector; using std :: ostream;

104

using std :: pair;

using std :: make_pair;

#include <CGAL/Delaunay_triangulation_2.h>

#include <CGAL/Triangulation_vertex_base_with_info_2.h>

const int NP = 2;

class TSection2{

friend ostream& operator<< (ostream &, const TSection2&);

public:

//TSection2(void);

TSection2(

const double,

const vector<Vector2>&, const vector<Vector2>&, const vector<double>&, const double _rho

);

const vector<Vector2>& GetX() const; const vector<Vector2>& GetP() const; const vector<double>& GetZ() const;

const TSection2& operator = (const TSection2&); bool operator == (const TSection2&) const; double Time()const{return t;}

void Formula();

const double GetRho()const{return Rho;}

private:

double t; vector<Vector2> X; vector<Vector2> P; vector<double> Z; const double Rho;

typedef CGAL::Triangulation_vertex_base_with_info_2<size_t, K> Vb;

typedef CGAL::Triangulation_data_structure_2<Vb>		Tds;
typedef CGAL::Delaunay_triangulation_2<K, Tds>		Triangulation;
typedef Triangulation::Point	Point;

Triangulation Triangulate(); void FindIntersectionAndMin();

void FindIntersectionAndMin(Triangulation &T);

};

//TSECTION2.CPP #include "TSection2.h"

#include "../Include/Functions.h" #include <math.h> #include<algorithm>

#include <map> #include <functional>

using std :: map; using std :: fill; using std :: swap;

using std :: mem_fun_ref;

ostream& operator<< (ostream &Out, const TSection2 &Section){ Out << "t = " << Section.Time() << '\n';

Out << "X = [\n";

Out << Section.GetX(); Out << "]\n";

Out << "P = [\n";

105

Out << Section.GetP(); Out << "]\n";

Out << "Z = [\n";

Out << Section.GetZ(); Out << "]\n";

return Out;

}

TSection2::TSection2( const double _t,

const vector<Vector2>& _X, const vector<Vector2>& _P, const vector<double>& _Z, const double _Rho

)

:t(_t), X(_X), P(_P), Z(_Z), Rho(_Rho)

{}

const vector<Vector2>& TSection2::GetX()const{ return X;

}

const vector<Vector2>& TSection2::GetP()const{ return P;

}

const vector<double>& TSection2::GetZ()const{ return Z;

}

const TSection2& TSection2 :: operator = (const TSection2 &Section){ t = Section.t;

X = Section.X;

P = Section.P;

Z = Section.Z; return *this;

}

bool TSection2 :: operator == (const TSection2 &Section) const { return X == Section.X

&&P == Section.P

&&Z == Section.Z;

}

TSection2::Triangulation TSection2 :: Triangulate()

{

std::vector< std::pair<Point,size_t> > points; for(size_t i=0; i < X.size();++i)

points.push_back( std::make_pair(CGAL::ORIGIN+X[i],i) ); Triangulation T;

//clock_t timer=clock();

//std::cout<<"Triangulation time: ";//38595 мсек - 1999 точек /// 1390000 мсек - 12500 точек

T.insert( points.begin(),points.end() ); //std::cout<<clock()-timer<<std::endl; return T;

}

void TSection2 :: FindIntersectionAndMin(TSection2::Triangulation &T)

{

vector<pair<size_t,size_t> > OldNewIndex;
Triangulation::Finite_vertices_iterator	it	= T.finite_vertices_begin(),

end =

T.finite_vertices_end();

size_t Size = T.number_of_vertices(); Triangulation::Finite_vertices_iterator* ArrayIterators

106

= new Triangulation::Finite_vertices_iterator[Size]; size_t i = 0;

for (; it != end; ++it){ ArrayIterators[i++] = it;

}

omp_lock_t Lock; omp_init_lock(&Lock);

#pragma omp parallel for num_threads(NP) for(int i = 0; i < Size; ++i)

{

assert(omp_get_num_threads() <= omp_get_num_procs()); assert(omp_get_num_threads() == NP);

Triangulation::Vertex_circulator vc = T.incident_vertices(ArrayIterators[i]); if(vc != 0){

Triangulation::Vertex_circulator done = vc; Point p = ArrayIterators[i]->point();

do {

if(T.is_infinite(vc)) continue;

Point p1 = vc->point();

if(vc->info() == ArrayIterators[i]->info())

{

std::cout<<"assert"<<std::endl;

}

CGAL :: Comparison_result Dist = CGAL::compare_squared_distance ( p, p1, Rho

); if (

Dist == CGAL :: SMALLER || Dist == CGAL :: EQUAL

){//если расстояние меньше RhoX if(Z[ArrayIterators[i]->info()] > Z[vc->info()])

{

omp_set_lock(&Lock); OldNewIndex.push_back(make_pair(ArrayIterators[i]-

>info(),vc->info()));

omp_unset_lock(&Lock);

}

}while(++vc != done);

}

omp_destroy_lock(&Lock); delete []ArrayIterators;

for(size_t i = 0; i < OldNewIndex.size(); ++i)

{

Z[OldNewIndex[i].first] = Z[OldNewIndex[i].second];

P[OldNewIndex[i].first] = P[OldNewIndex[i].second];

}

void TSection2 :: FindIntersectionAndMin()

{

vector<double> vMin = Z; for(size_t i = 0; i < X.size(); ++i)

{

for(size_t j = i+1; j < X.size(); ++j)

{

if(X[i] == X[j])

{

if(vMin[i] > Z[j])

{

vMin[i] = Z[j];

107

}

if(vMin[j] > Z[i])

{

vMin[j] = Z[i];

}

Z = vMin;

}

void TSection2 :: Formula()

{

//FindIntersectionAndMin(); Triangulation temp = Triangulate(); FindIntersectionAndMin(temp);

}

//TVALUEFUNCTION.H #pragma once

#include <fstream> #include <iostream> #include <vector> #include <algorithm> #include "TSection2.h" #include "Tree2.h" #include "../Include/Tree.h" using std :: string;

using std :: ofstream; using std :: ifstream; using std :: ostream; using std :: istream; using std :: endl;

class TValueFunction2

{

public:

TValueFunction2

(

const TTree &times,

const vector<TTree2> &Trajectories, const vector<TTree> &values

)

:Times(times), Trajectories(Trajectories), Values(values)

{}

void Instant2Files(const string &FileName); bool operator == (const TValueFunction2 &V1);

friend ostream& operator << (ostream &Out, const TValueFunction2 &VF); vector<Vector2> GetSectionPoints(const size_t SectionNum) const;

private:

TTree Times; vector<TTree2> Trajectories; vector<TTree> Values;

};

//TVALUEFUNCTION.CPP #include "TValueFunction2.h" #include "../Include/Functions.h" #include "Problem2.h"

#include <assert.h> using std :: ofstream; using std :: string;

using std :: ostringstream;

108

bool TValueFunction2::operator == (const TValueFunction2 &V1)

{

return V1.Trajectories == this->Trajectories;

}

ostream& operator << (ostream &Out, const TValueFunction2 &VF)

{

Out << VF.Trajectories; return Out;

}

void TValueFunction2::Instant2Files(const string &DirName){ assert(!Times.Empty());

unsigned int Size = Times.Size(); for(unsigned int t = 0; t < Size; ++t){

ostringstream FileName;

FileName << DirName+"/t_" << t << ".txt"; ofstream Out(FileName.str().c_str());

Out << 3 << " " << Trajectories.size() << endl; for(unsigned int k = 0; k < Trajectories.size(); ++k){

Out << Trajectories[k][t].x() << " "; Out << Trajectories[k][t].y() << " "; Out << Values[k][t] << endl;

}

Out.close();

}

vector<Vector2> TValueFunction2::GetSectionPoints(const size_t SectionNum) const

{

vector<Vector2> Result;

for(size_t i = 0; i < Trajectories.size(); ++i)

{

Result.push_back(Trajectories[i].GetBranche(SectionNum));

}

return Result;

}

109

<<< < Предыдущая 1 2 3 4 5 6 7 89 / 99

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

#
22.02.20152.64 Mб1155Diplom.doc
#
22.02.2015296.45 Кб14diploma_method.doc
#
22.02.2015358.91 Кб18Diploma_Thesis-P_Adamec.doc
#
22.02.20153 Mб88diplom_pz (для примера оформления).doc
#
21.09.20198.62 Mб13Diplom_PZ.doc
#
23.02.20157.89 Mб17Diplom_SamoylovAL.pdf
#
13.03.20164.14 Mб7DiscDogsRockEnglish.pdf
#
13.03.2016383.48 Кб80dist_8.pdf
#
27.09.2019122.63 Кб11DKB_bilety.docx
#
02.05.2019304.19 Кб3dlya_061.docx
#
24.09.2019913.41 Кб8Dlya_studentov-LEKTsII.doc