Диссертация (1155087), страница 10
Текст из файла (страница 10)
и, проинтегрировав полученное равенство по переменной от 0 до 1 придем к выражению:∑︁∑︁ ∑︁= () − () − ()2 .=1=1=1(А.5)106Окончательно получим:⎧∑︀∑︀⎪0⎪0⎪=()− 0 ()⎨ ⎪⎪⎪⎩ ==1∑︀ ()−=1=1∑︀ ()−=1∑︀(А.6)2 () .=1После интегрирования системы (А.6) восстанавливались неизвсетные функции (, ).Далее приведен класс, разработанные на языке C++, позволяющий находить численное решение распределенной математической модели квазивидов:определение класса#include <list>//для передачи начальной функкции u(x,0)typedef double (*myfunc)(double /*x*/, int /*i*/, int /*n*/);//Далее везде считается, что для матриц: C(i , j): i - номер неизвестной функции,// j - номер приближения//Размерность: С(n,m); (n - уранений, m - приближений)class SolVer{public://получить решениеvoid get_solution(myfunc/*u_initial - начальное значение u*/,Matrix<double> Q/*матрица вероятностей */,Matrix<double> /*a - матрица коэффициентов*/,Matrix<double> /*D - матрица диффузии*/,char* /*name_file - файл, в который записать решение*/);//инициализировать значение шагов по времени и координате//m - размерность апроксимацииSolVer(int /*n*/, int /*m*/, double /*ht*/, double /*xa*/,double /*xb*/, double /*ta*/, double /*tb*/);private://Получить начальные условия для C107Matrix<double> get_initial_value(myfunc/*u_initial - начальное значение u*/,int /*n*/ , int /*m*/);//записать начальное значение uvoid write_initial_value(char* /*name_file*/);//получить фитнессdouble get_fitness(Matrix<double> /*a*/);//сумма по k = 1,...,n:a(i , k) * C(k , j);double get_summC_with_koeff(Matrix<double> /*a*/, Matrix<double> /*C*/,int /*i*/, int /*j*/);//индикатор: записывать или нет значениеbool indicator_of_write(double /*t*/, double /*x*/);double get_dissipation(Matrix<double> /*D*/, Matrix<double> /*C*/,int /*i*/, int /*j*/);private://коэффициенты разложения (u(i) = \sum_{j=1}^n c(i,j)*cos((j-1)*pi*x)Matrix<double> C;Matrix<double> u; //решение системыprivate:double ht; //шаг по временияdouble hx; //шаг по координатеdouble xa; //левая граница решения по координатной переменнойdouble xb; //правая граница решения по координатной переменнойdouble ta; //левая граница решения по временной переменнойdouble tb; //правая граница решения по временной переменнойint n;//количество уравнений в частных производныхint m;//размерность аппроксимацииdouble delta_xb; //решать пока x < xb + delta_xb};108реализация класса#include "stdafx.h"#include "SolVer.h"#include <fstream>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);this->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++)109{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;for(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())110output << u(i) << ’\t’;}output << u(C.readn()) << std::endl;u = 0.0;}}void SolVer::get_solution(myfunc u_initial, Matrix<double> q, Matrix<double> a,Matrix<double> D, char* name_file){C = get_initial_value(u_initial, n, m);Matrix<double> qa = q * a;write_initial_value(name_file);std::ofstream output(name_file, std::ios_base::app); //дозаписывать файлfor(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++){C(i, 1) += ht * (get_summC_with_koeff(qa, C, i, 1) get_fitness(qa) * C(i,1));u(i) = u(i) + C(i , 1);//readm - размерность апроксимацииfor(int l = 2; l != C.readm() + 1; l++){C(i , l) += ht * (get_summC_with_koeff(qa, C, i, l) get_fitness(a) * C(i , l) -111(l * pi * l * pi) * get_dissipation(D, C, i, l));u(i) = u(i) + C(i , l) * cos((l - 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_summC_with_koeff(Matrix<double> qa, Matrix<double> C,int i, int l){double S = 0.0;Matrix<double> C0(C.readn());for(int j = 1; j != C0.readn() + 1; j++)C0(j) = C(j, l);Matrix<double> qc = qa * C0;return qc(i);}bool SolVer::indicator_of_write(double t, double x){//if(t < 1 || (t > 1 && (int(10*t)%10 ==0)))// if(x == 0.5)//if(int(10*t)%10 ==0 && t > 1)112return true;return false;}double SolVer::get_fitness(Matrix<double> qa){double fs = 0.0;Matrix<double> C0(C.readn());for(int j = 1; j != C0.readn() + 1; j++)C0(j) = C(j, 1);Matrix<double> qf = qa * C0;for(int i = 1; i != qa.readn() + 1; i++)fs += qf(i);//C(i , 1) * qa(i, i);return fs;}double SolVer::get_dissipation(Matrix<double> D, Matrix<double> C, int i, int l){double S = 0.0;for(int j = 1; j != D.readn() + 1; j++){S += D(i,j) * C(j, l);}return S;}113Приложение БЧисленное решение математической моделидвойного гиперциклаУравнение двойного гипрцикла:˙ = ( −1 −1 −2 − ),где=∑︁ = 1, ,(Б.1) −1 −1 −2=1Выполняются следующие ограничения: ∈ ; = { | = 1, ⇒ ≥ 0;∑︀=1 = 1}(Б.2)При численном решении системы (Б.1) необходимомо, что бы выполненялось ограничение (Б.2) для любого значения .Теорема Б.1.При аппроксимации методом Рунге - Кутта произвольного порядка относительно системы(Б.1) ограничения (Б.2) выполняются автоматически.ДоказательствоСогласно методу Рунге - Кутта для -го уравненияимеем+1=+ℎ∑︁ ,(Б.3)=1где , - весовые коэффициенты.Просуммируем левую и правую часть (Б.3) по :=∑︁=1∑︁∑︁+ℎ =1=1(Б.4)114Первое слагаемое суммы∑︀=1 = 1.
Рассмотрим теперь второе слагаемоесуммы (Б.3) и покажем, что оно равно 0: ∑︁∑︁∑︁∑︁ = ℎℎ=1(Б.5)=1 =1=1При = 1 сумма (Б.5) равнаℎ1∑︁1 = ℎ1∑︁=1( , ) = ℎ1=1∑︁(︃(︃−1 −2 −=1∑︁)︃)︃ −1 −2=0=1(Здесь ( , ) - правая часть уравнения 1)При = 2:∑︀∑︀2 = =1 ( + 2 ℎ, + 21 ℎ1 ) ==1(︁)︁)︁ (︁∑︀=( + 21 ℎ1 ) (−1 + 21 ℎ1 )(−2 + 21 ℎ1 ) − ,=1где=∑︁По доказанномупреобразуется в∑︁( + 21 ℎ1 )(−1 + 21 ℎ1 )(−2 + 21 ℎ1 )=1∑︀( + 21 ℎ1 ) = 1, следовательно, наше выражение=1( + 21 ℎ1 )(−1 + 21 ℎ1 )(−2 + 21 ℎ1 ) − = 0=1(по определению ).
Допустим, что мы доказали выполнение условий 2для − 1-го приближения. Докажем для -го:∑︀ ==1∑︀( + ℎ, + 1 ℎ1 + ... + ,−1 ℎ−1 ) ==1= (( + 1 ℎ1 + ... + ,−1 ℎ−1 )·=1·(((−1 + 1 ℎ1 + ... + ,−1 ℎ−1 )((−2∑︀∑︀+ 1 ℎ1 + ... + ,−1 ℎ−1 ) − )),=1 (( + 1 ℎ1 + ... + ,−1 ℎ−1 )(−1 + 1 ℎ1 + ...∑︀,−1 ℎ−1 )(−2 + 1 ℎ1 + ... + ,−1 ℎ−1 )). Но из того, что (=1где =++1151 ℎ1 + ... + ,−1 ℎ−1 ) = 1 ⇒∑︀ = 0=1Это и означает, что применение метода Рунге-Кутта здесь корректно.Численное решение распределенной математической моделидвойного гиперциклаДано уравнение:(︃= −1 −1 −2 −1 R∑︀1 R∑︀)︃ −1 −1 −2 =1 0 (, ) = 1;=1 0 (, 0) = (); (0, )= (1, )+ 2 (,)2 ,(Б.6)= 0,() ≥ 0.Решение (, ) будем искать в виде (, ) =∑︁(Б.7) () cos()=0В этом случае получаем:∑︀=0= ()∑︀cos() = () cos()( −1=0=0+ ·∑︀∑︀−1 () cos())∑︀−2 () cos()) − +=0(−) () cos(),=0где =∑︀(Б.8) · −1 −2 ; − Выражение от .Проинтегрируем последнее выражение по от 0 до 1, получим системууравнений для определения неизвестных функций ().
Получившаясясистема достаточно громоздкая. Заметим, что для решения этой системынеобходимо вычисление интгерала от произведения косинусов. Для этогобыл разработан класс на языке C++, позволяющий вычислять интегралот произведения cos(1 ) · ... · cos( ) при целых и ≤ 4:116Определение класса:class Numerical_calc_integral{public:Numerical_calc_integral();//вычисление итнеграла от cos(k1 * pi * x), по переменной x на интервале// от 0 до 1double Numerical_calc_cos(int /*k1*/);//вычисление итнеграла от cos(k1 * pi * x) * cos(k2 * pi * x),//по переменной x на интервале от 0 до 1double Numerical_calc_cos(int /*k1*/, int /*k2*/);//... cos(k1 * pi * x) * cos(k2 * pi * x) * cos(k3 * pi * x)double Numerical_calc_cos(int /*k1*/, int /*k2*/, int /*k3*/);//...