math-library

1

2
#include <assert.h>
3
#include <iostream>
4
#include <fstream>
5
#include <math.h>
6

7
#include "Constants.h"
8
#include "Engine.h"
9

10
#define L(x1,x2)	(fabs(x1-x2))
11

12
using namespace std;
13

14
void Progonka(int N, CMatrix &A,CMatrix &B,CMatrix &C,CMatrix &F, double m1,double n1,double m2,double n2, CMatrix &Y)
15
{
16
	CMatrix Alpha(N+1), Beta(N+1);
17

18
	Alpha(1) = m1; Beta(1) = n1;
19
	for(int i=1; i<N; i++)
20
	{
21
		Alpha(i+1) = B(i)/(C(i) - A(i)*Alpha(i));
22
		Beta(i+1) = (A(i)*Beta(i) + F(i))/(C(i) - A(i)*Alpha(i));
23
	}
24

25
	Y(N) = (n2 + m2*Beta(N))/(1 - m2*Alpha(N));
26
	for(int i=N-1; i>=0; i--) Y(i) = Alpha(i+1)*Y(i+1) + Beta(i+1);
27
}
28

29
double Simpson(CInterval &AB, double (*getF)(double))
30
{
31
	if(AB.N() % 2) AB.ReBorn(AB.X1(), AB.X2(), AB.N()+1);	//N должно быть четным
32

33
	double I = 0;
34

35
	for(int i=1; i<AB.N(); i+=2) I += getF(AB.X(i));
36
	I *= 2;
37
	for(int i=2; i<=AB.N()-2; i+=2) I += getF(AB.X(i));
38
	I = AB.H()/3.*(
39
		getF(AB.X1()) + I*2 + getF(AB.X2())
40
		);
41

42
	return I;
43
}
44

45
void mSOLVE(CMatrix &a, CMatrix &b, CMatrix &X)
46
{
47
	int Size = a.GetM();
48
	assert(a.GetN() == Size && b.GetM() == Size);
49

50
	CMatrix	A(a), B(b);		//Создаем копии матриц, т.к. над ними проводим действия
51

52
	//Приведение матрицы к треугольному виду
53
	for(int Curr=0; Curr<Size; Curr++)
54
	{
55
		int MaxNm = Curr;
56

57
		//Поиск максимального элемента в столбце Curr
58
		for(int i=Curr+1; i<Size; i++)
59
			if(A(i,Curr)>A(i,MaxNm)) MaxNm = i;
60

61
		//Перестановка строк
62
		if(MaxNm!=Curr)
63
		{
64
			A.SwapLines(Curr,MaxNm); B.SwapLines(Curr,MaxNm);
65
		}
66

67
		//Деление на коэффициент при первом члене
68
		B(Curr) = B(Curr)/A(Curr,Curr);
69
		for(int i=Size-1; i>=Curr; i--)	A(Curr,i) = A(Curr,i)/A(Curr,Curr);
70

71
		//Вычитание строк
72
		for(int line=Curr+1; line<Size; line++)
73
		{
74
			B(line) = B(line) - B(Curr)*A(line,Curr);
75
			for(int i=Size-1; i>=Curr; i--)
76
				A(line,i) = A(line,i) - A(Curr,i)*A(line,Curr);
77
		}
78
	}
79

80
	//Нахождение решения системы по треугольной матрице
81
	double d;
82
	for(int i=Size-1; i>=0; i--)
83
	{
84
		d = B(i);
85
		for(int j=i+1; j<Size; j++) d -= A(i,j)*X(j);
86
		X(i) = d/A(i,i);
87
	}
88
}
89

90
void mINV(CMatrix &A, CMatrix &invA)
91
{
92
	assert(A.GetN() == A.GetM());
93

94
	int Size = A.GetM();
95
	invA.SetSize(Size,Size);
96

97
	CMatrix	E(Size),T(Size);
98
	for(int i=0; i<Size; i++)
99
	{
100
		E = 0; E(i) = 1;
101
		mSOLVE(A,E,T);
102

103
		for(int j=0; j<Size; j++) invA(j,i) = T(j);
104
	}
105
}
106

107
double mDET(CMatrix &A)
108
{
109
	int Size = A.GetM();
110
	assert(A.GetN() == Size);
111

112
	int NSw = 0; //Количество перестановок строк
113

114
	//Приведение матрицы к треугольному виду
115
	for(int Curr=0; Curr<Size; Curr++)
116
	{
117
		int MaxNm = Curr;
118

119
		//Поиск максимального элемента в столбце Curr
120
		for(int i=Curr+1; i<Size; i++) if(A(i,Curr)>A(i,MaxNm)) MaxNm = i;
121

122
		//Перестановка строк
123
		if(MaxNm!=Curr)	{ A.SwapLines(Curr,MaxNm); NSw++; }
124

125
		//Деление на коэффициент при первом члене
126
//		for(i=Size-1; i>=Curr; i--)	A(Curr,i) = A(Curr,i)/A(Curr,Curr);			
127

128
		//Вычитание строк
129
		for(int line=Curr+1; line<Size; line++)
130
		for(int i=Size-1; i>=Curr; i--)	A(line,i) = A(line,i) - A(Curr,i)/A(Curr,Curr)*A(line,Curr);
131
	}
132

133
	double DET = (NSw % 2) ? -1:1;
134
	for(int i=0; i<Size; i++) DET *= A(i,i);
135

136
	return DET;
137
}
138

139
double FindMinMax(CInterval &AB, double(*f)(double), int find_max)
140
{
141
	double x[4];
142
	x[0] = AB.X1(); x[2] = AB.X2();
143
	x[1] = x[0] + 2/(3+sqrt(5))*(x[2]-x[0]);
144

145
	while(x[2]-x[0]>AB.H())
146
	{
147
		x[3] = x[0]+x[2]-x[1];	
148
		int	m_i = 0, f_i = 0;
149
	
150
		//поиск минимальной из 4 точек
151
		for(int i=1; i<4; i++) if(find_max*(f(x[i])-f(x[m_i])) > 0) m_i = i;
152

153
		//поиск наиболее удаленной от минимума точки
154
		for(int i=1; i<4; i++) if(L(x[m_i],x[i]) > L(x[m_i],x[f_i])) f_i = i;
155

156
		//перенумеровываем точки, чтобы наиболее удаленная стала 4-й
157
		if(f_i!=3) x[f_i] = x[3];
158

159
		//упорядочение
160
		double d;
161

162
		for(m_i=0; m_i<=1; m_i++)
163
		{
164
			int min = m_i;
165
			for(int i=1; i<=2; i++) if(x[i]<x[m_i]) min = i;
166

167
			d = x[min];
168
			x[min] = x[m_i];
169
			x[m_i] = d;
170
		}
171
	}
172
	return 0.5*(x[2]+x[0]);
173
}
174

175
void SaveFunc(CInterval &AB, double (*func)(int), char *fname)
176
{
177
	ofstream f(fname);
178
	for(int i=0; i<=AB.N(); i++) f<< AB.X(i) <<" "<< func(i) <<"\n";
179
	f.close();
180
}
181

182
void SaveFunc(CInterval &AB, double (*func)(double), char *fname)
183
{
184
	ofstream f(fname);
185
	for(int i=0; i<=AB.N(); i++) f<< AB.X(i) <<" "<< func(AB.X(i)) <<"\n";
186
	f.close();
187
}
188

189
void SaveFunc(CInterval &t, double(*x)(double), double (*func)(double), char *fname)
190
{
191
	ofstream f(fname);
192
	for(int i=0; i<=t.N(); i++) f<< x( t.X(i) )<<" "<< func( t.X(i) ) <<"\n";
193
	f.close();
194
}
195

196
void SaveCorrelationFunc(CMatrix &A, char *fname, double precision)
197
{
198
	assert(0<precision && precision<=1 && A.GetM()>=1);
199
	
200
	double Min = A.Min(), Max = A.Max();
201
	CInterval AB(Min, Max, (Max-Min)*precision);
202
	CMatrix P(AB.N(),2);			//Частоты
203
	P = 0;
204
	
205
	for(int i=0; i<A.GetM(); i++)
206
	{
207
		double f_i = A(i);
208
		int num_of_f_i = AB.i(f_i);
209
		
210
		if(num_of_f_i != AB.N())	//CInterval включает обе границы диапазона 0..N; CMatrix - только 0
211
		{
212
			P(num_of_f_i,0) = A(i);
213
			P(num_of_f_i,1) += 1;
214
		}
215
	}
216
	
217
	//Нормировка
218
	//
219
	double max = P(0,1);
220
	for(int i=1; i<P.GetM(); i++) if(P(i,1) > max) max = P(i,1);
221
	for(int i=0; i<P.GetM(); i++) P(i,1) /= max;
222
	
223
	P.Save(fname);
224
}
225

226
void SaveCorrelationFunc(CInterval &AB, double (*func)(double), char *fname, double precision)
227
{
228
	CMatrix A(AB.N());
229
	for(int i=0; i<AB.N(); i++) A(i) = func(AB.X(i));
230
	
231
	SaveCorrelationFunc(A,fname,precision);
232
}
233

234
CMatrix::CMatrix(int m, int n)
235
{
236
	assert(m>=1 && n>=0);
237
	pData = new double[(M = m)*(N = n)];
238
	assert(pData);
239
}
240

241
void CMatrix::Save(char *fname, bool orig_view)
242
{
243
	ofstream f(fname);
244

245
	if(orig_view)
246
		for(int i=0; i<M; i++)
247
		{
248
			for(int j=0; j<N; j++)	f<< (*this)(i,j) <<" ";
249
			f<<"\n";
250
		}
251
	else
252
		for(int j=0; j<N; j++)
253
		{
254
			for(int i=0; i<M; i++)	f<< (*this)(i,j) <<" ";
255
			f<<"\n";
256
		}
257

258
	f.close();
259
}
260

261
void CMatrix::SwapLines(int m1, int m2)
262
{
263
	assert(m1<M && m2<M && m1>=0 && m2>=0);
264
	
265
	double el;
266
	for(int i=0; i<N; i++)
267
	{
268
		el = *(pData + m1*N + i);
269
		*(pData + m1*N + i) = *(pData + m2*N + i);
270
		*(pData + m2*N + i) = el;
271
	}
272
}
273

274
void CMatrix::SwapCols(int n1, int n2)
275
{
276
	assert(n1<N && n2<N && n1>=0 && n2>=0);
277

278
	double el;
279
	for(int i=0; i<N; i++)
280
	{
281
		el = *(pData + i*N + n1);
282
		*(pData + i*N + n1) = *(pData + i*N + n2);
283
		*(pData + i*N + n2) = el;
284
	}
285
}
286

287
CMatrix &CMatrix::operator=(CMatrix &right)
288
{
289
	if(&right != this)
290
	{
291
		SetSize(right.M,right.N);
292

293
		for(int i=0; i<M; i++)
294
		for(int j=0; j<N; j++) *(pData + i*N + j) = *(right.pData + i*N + j);
295
	}
296
	return *this;
297
}
298

299
CMatrix &CMatrix::operator=(double d)
300
{
301
	for(int i=0; i<M; i++)
302
	for(int j=0; j<N; j++) *(pData + i*N + j) = d;
303
		
304
	return *this;
305
}
306

307
CMatrix &CMatrix::operator-()
308
{
309
	for(int i=0; i<M; i++)
310
	for(int j=0; j<N; j++) *(pData + i*N + j) = -(*(pData + i*N + j));
311
		
312
	return *this;
313
}
314

315
bool CMatrix::operator==(CMatrix &right)
316
{
317
	assert(M==right.M && N==right.N);
318

319
	for(int i=0; i<M; i++)
320
	for(int j=0; j<N; j++) if(*(pData + i*N + j) != *(right.pData + i*N + j)) return false;
321
	
322
	return true;
323
}
324

325
double Round(double d)
326
{
327
	double r = ceil(d);		//Использование (int) недопустимо!
328
	if(fabs(d-r) > 0.5) return r += sign(d); else return r;
329
}
330

331
CMatrix::~CMatrix()
332
{
333
	if(!pData) delete [] pData;
334
}
335

336
void CMatrix::SetSize(int m, int n)		//!Память не освобождается, если L-массив УЖЕ имеет нужные размеры
337
{
338
	assert(m>=1 && n>=0);
339

340
	if(M!=m || N!=n)
341
	{
342
		if(!pData) delete [] pData;
343
		pData = new double[(M = m)*(N = n)];
344
		assert(pData);
345
	}
346
}
347

348
double CMatrix::Max()
349
{
350
	double Max = (*this)(0,0);
351
	
352
	for(int i=0; i<M; i++)
353
	for(int j=0; j<N; j++) if((*this)(i,j)>Max) Max = (*this)(i,j);
354

355
	return Max;
356
}
357

358
double CMatrix::Min()
359
{
360
	double Min = (*this)(0,0);
361

362
	for(int i=0; i<M; i++)
363
	for(int j=0; j<N; j++) if((*this)(i,j)<Min) Min = (*this)(i,j);
364

365
	return Min;
366
}
367

368
CMatrix::CMatrix(CMatrix &M)
369
{
370
	*this = M;
371
}
372

373
void empty_func(CLattice &Lattice) {}
374

375
void MMD(CLattice &Lattice, double dt, int StepsCount, bool sEnergy, void (*func)(CLattice &))
376
{
377
	int i,k, N = Lattice.N;
378
	CMatrix F1(N,3),F2(N,3);
379
	double t_2 = dt/2.;
380
	bool flag;
381
	
382
	Lattice.CalcForces(F1);
383
	
384
	double E0 = Lattice.getE_k() + Lattice.getU();	//Потенциальная энергия известна только после вызова CalcForces()
385

386
	for(int step=0; step<StepsCount; step++)
387
	{
388
		for(k=0; k<3; k++)
389
		for(i=0; i<N; i++) Lattice.XYZ(i,k) += dt*(Lattice.VXYZ(i,k) + t_2*F1(i,k));
390
		
391
		flag = CheckValue(Lattice.XYZ,Lattice.L);	//Вылетевшие частицы впускаем через противоположную стенку
392
		if(!flag)									//Проверка, не пролетает ли частица за dt больше L
393
		{
394
			cout<< "\nstep = "<< step <<"\n! dt is large !\n" <<endl; exit(1);
395
		}
396
		
397
		Lattice.CalcForces(F2);
398
			
399
		for(k=0; k<3; k++)
400
		for(i=0; i<N; i++) Lattice.VXYZ(i,k) += t_2*(F1(i,k) + F2(i,k));
401
		F1 = F2;
402

403
		if(sEnergy)
404
		{
405
			double H = 0;
406
			for(k=0; k<3; k++)
407
			for(i=0; i<N; i++) H += pow(Lattice.VXYZ(i,k),2) + pow(F2(i,k),2);
408

409
			double dE_H = (Lattice.getE_k()+Lattice.getU() - E0)/H;
410

411
			for(k=0; k<3; k++)
412
			for(i=0; i<N; i++)
413
			{
414
				Lattice.XYZ(i,k) += F2(i,k)*dE_H;
415
				Lattice.VXYZ(i,k) -= Lattice.VXYZ(i,k)*dE_H;
416
			}
417
		} //if
418
		
419
		func(Lattice);
420
	}
421
}