Диссертация (1155087), страница 11
Текст из файла (страница 11)
cos(k1 * pi * x) * cos(k2 * pi * x) * cos(k3*pi*x) * cos(k4*pi*x)double Numerical_calc_cos(int /*k1*/, int /*k2*/, int /*k3*/, int /*k4*/);};Реализация класса:Numerical_calc_integral::Numerical_calc_integral(){}double Numerical_calc_integral::Numerical_calc_cos(int k){if(k == 0)return 1.0;return 0.0;}117double Numerical_calc_integral::Numerical_calc_cos(int k1, int k2){if(k1 == 0)return Numerical_calc_cos(k2);if(k2 == 0)return Numerical_calc_cos(k1);if(k1 == k2)return 0.5;return 0.0;}double Numerical_calc_integral::Numerical_calc_cos(int k1, int k2, int k3){if(k1 == 0)return Numerical_calc_cos(k2, k3);if(k2 == 0)return Numerical_calc_cos(k1, k3);if(k3 == 0)return Numerical_calc_cos(k1, k2);if((k1 == k2 + k3) || (k2 == k1 + k3) || (k3 == k1 + k2))return 0.25;return 0.0;}double Numerical_calc_integral::Numerical_calc_cos(int k1, int k2, int k3, int k4){if(k1 == 0)return Numerical_calc_cos(k2, k3, k4);if(k2 == 0)return Numerical_calc_cos(k1, k3, k4);if(k3 == 0)118return Numerical_calc_cos(k1, k2, k4);if(k4 == 0)return Numerical_calc_cos(k1, k2, k3);if((k1 == k2 + k3 + k4) || (k2 == k1 + k3 + k4) ||(k3 == k1 + k2 + k4) || (k4 == k1 + k2 + k3))return 0.125;if((((k1 == k2)&& (k3 == k4)) && (k1 != k3)) || (((k1 == k3)&& (k2 == k4))&& (k2 != k3)) || (((k1 == k4)&& (k2 == k3)) && (k1 != k2)))return 0.25;if((k1 == k2) && (k2 == k3) && (k3 == k4))return 0.375;if((k1 + k2 == k3 + k4) || (k1 + k3 == k2 + k4) || (k1 + k4 == k2 + k3))return 0.125;return 0.0;}Далее приведен класс, разработанный на языке C++, позволяющий находить численное решение рапспределенной математической модели двойного гиперцикла:определение класса#include "stdafx.h"#include "CNumerical_calc_integral.h"#include <iostream>#include <list>//для передачи начальной функкции u(x,0)typedef double (*myfunc)(double /*x*/, int /*i*/, int /*n*/);//Далее везде считается, что для матриц: C(i , j):// i - номер неизвестной функции, j - номер приближения//Размерность: С(n,m); (n - уранений, m - приближений)119class SolVer{public://получить решениеvoid get_solution(myfunc/*u_initial - начальное значение u*/, Matrix<double>/*D - матрица диффузии*/,char* /*input_name - файл, в который записать решение*/);//инициализировать значение шагов по времени и координате//m - размерность апроксимацииSolVer(int /*n*/, int /*m*/, double /*ht*/, double /*xa*/,double /*xb*/, double /*ta*/, double /*tb*/);private://получить фитнессdouble get_fitness();//первый интеграл в суммеdouble get_integral1(int /*i - номер уравнения*/,int /*j - домножаем на cos(j*pi*x)*/);//Получить начальные условия для CMatrix<double> get_initial_value(myfunc/*u_initial - начальноезначение u*/,int /*n*/ , int /*m*/);void write_initial_value(char* ); //записать начальное значение uvoid write_new_value(double /*t*/); //записать очередное значение в файл//индикатор: записывать или нет значениеbool indicator_of_write(double /*t*/, double /*x*/);private:char* name_file;120private:Matrix<double> C; //коэффициенты разложенияMatrix<double> u; //решение системыNumerical_calc_integral Integral; //для вычисления интегралаprivate:double ht; //шаг по временияdouble hx; //шаг по координатеdouble xa; //левая граница решения по координатной переменнойdouble xb; //правая граница решения по координатной переменнойdouble ta; //левая граница решения по временной переменнойdouble tb; //правая граница решения по временной переменнойint n;//количество уравнений в частных производныхint m;//размерность аппроксимацииdouble delta_xb; //решать пока x < xb + delta_xb//===================private:Matrix<double> k1(Matrix<double> );Matrix<double> k2(Matrix<double> );Matrix<double> k3(Matrix<double> );};реализация класса#include "stdafx.h"#include "SolVer.h"#include <fstream>#include <iomanip>SolVer::SolVer(int n, int m, double ht, double xa,double xb, double ta, double tb): C(Matrix<double> (n , m)), u(Matrix<double> (n)){this->ht = ht;this->hx = (xb - xa)/(m - 1);121this->xa = xa;this->xb = xb;this->ta = ta;this->tb = tb;this->n = n;this->m = m;delta_xb = xb / 10000.0;u = 0.0;}//получить начальные значения для неизвестных функций C(i,j)//n - количество уравнений в частных производных//m - размерность аппроксимацииMatrix<double> SolVer::get_initial_value(myfunc u_initial, int n , int m){//матрица коэффициентов СЛАУ для отыскания неизвестных функциийMatrix<double> A(m , m);//вектор правых частей СЛАУ для отыскания неизвестных функциийMatrix<double> B(m);for(int j = 1; j != n + 1; j++){int i = 1;for(double x = xa; x < xb + delta_xb; x += hx){for(int k = 1; k != m + 1; k++)A(k, i) = cos((k - 1) * pi * x);B(i) = (*u_initial)(x,j,n);i++;}//std::cout << B << std::endl;Matrix<double> C_j = MathLib::inv(A) * B;122for(int ii = 1; ii != m + 1; ii++)C(j , ii) = C_j(ii);}return C;}void SolVer::write_initial_value(char *name_file){std::ofstream output(name_file);for(double x = xa; x < xb + delta_xb; x += hx) //просчитать для всех x{output << 0 << ’\t’ << x << ’\t’;for(int i = 1; i != C.readn() + 1; i++){for(int j = 1; j != C.readm() + 1; j++)u(i) = u(i) + C(i , j) * cos((j - 1) * pi * x);if(i != C.readn())output << u(i) << ’\t’;}output << u(C.readn()) << std::endl;u = 0.0;}}void SolVer::get_solution(myfunc u_initial, Matrix<double> D, char* name_file){C = get_initial_value(u_initial, n, m);std::ofstream output(name_file, std::ios_base::app); //дозаписывать файл123for(double t = ta; t < tb; t+= ht) //просчитать для всех tfor(double x = xa; x < xb + delta_xb; x += hx) //просчитать для всех x{if(indicator_of_write(t, x))output << t + ht << ’\t’ << x << ’\t’;//readn - количество неизвестных функцийfor(int i = 1; i != C.readn() + 1; i++){//readm - размерность апроксимацииfor(int j = 1; j != C.readm() + 1; j++){C(i , j) = C(i , j) + ht* ((get_integral1(i, j - 1) -C(i , j) * Integral.Numerical_calc_cos(j - 1 , j - 1) *get_fitness() + D(i) * C(i,j) * (-(j-1)*pi*(j-1)*i) *Integral.Numerical_calc_cos(j - 1,j - 1)) /Integral.Numerical_calc_cos(j -1, j - 1));u(i) = u(i) + C(i , j) * cos((j - 1) * pi * x);}if(i != C.readn())if(indicator_of_write(t, x))output << u(i) << ’\t’;}if(indicator_of_write(t, x))output << u(C.readn()) << std::endl;u = 0.0;}}double SolVer::get_integral1(int i, int j){double S = 0.0;124int i_1 = (i==1)?C.readn():i-1;int i_2 = (i==2)?C.readn():i_1-1;for(int k2 = 1; k2 != C.readm() + 1; k2++)for(int k3 = 1; k3 != C.readm() + 1; k3++)for(int k4 = 1; k4 != C.readm() + 1; k4++){S += C(i, k2) * C(i_1, k3) * C(i_2, k4) *Integral.Numerical_calc_cos(k2 - 1, k3 - 1, k4 - 1, j);}return S;}double SolVer::get_fitness(){double fitness = 0.0;Numerical_calc_integral Integral;for(int i = 1; i != C.readn() + 1; i++){int i_1 = (i==1)?C.readn():i-1;int i_2 = (i==2)?C.readn():i_1-1;for(int k1 = 1; k1 != C.readm() + 1; k1++)for(int k2 = 1; k2 != C.readm() + 1; k2++)for(int k3 = 1; k3 != C.readm() + 1; k3++){fitness += C(i,k1) * C(i_1,k2) * C(i_2,k3) *Integral.Numerical_calc_cos(k1 - 1, k2 - 1, k3 - 1);}}return fitness;}125bool SolVer::indicator_of_write(double t, double x){//if(x ==0.5)if(int(10*t)%10 ==0 && t > -10)return true;return false;}.