Решение систем нелинейных уравнений https://www.mapleprimes.com/users/one%20man/posts?page=1
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

sysdiffeqn.cpp 9.6KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253
  1. #include "sysdiffeqn.h"
  2. #include <regex>
  3. //объявление статичных переменных
  4. std::vector<int> sys_diff_equation::vp_num ={};
  5. std::vector<int> sys_diff_equation::circ_num ={};
  6. state_type sys_diff_equation::circ_coord ={};
  7. unsigned sys_diff_equation::dim_line=2;
  8. std::vector<GiNaC::symbol> sys_diff_equation::var_name ={};
  9. sys_diff_equation::sys_diff_equation(sys_osn_dan &sod_): sod(&sod_)
  10. {
  11. std::string str, str_sm, dimension;
  12. std::stringstream _isstemp, mystreamin;
  13. std::vector<GiNaC::symbol> uprav_par_ex;
  14. // for (std::string str; std::getline(mystreamin_, str);)
  15. // if (!str.empty())
  16. // mystreamin << str << '\n';
  17. N=sod->get_par_diff().N;
  18. dt=sod->get_par_diff().dt;
  19. dim_line=sod->get_par_diff().dln;
  20. str_sm=sod->get_par_diff().str_sm;
  21. // mystreamin>>str_sm>>dt>>N>>dimension;
  22. // dim_line=dimension[0]-'0';
  23. // std::getline(mystreamin, str);
  24. // std::getline(mystreamin, str);
  25. x=sod->get_par_diff().x0;
  26. uprav_par_naz=sod->get_par_diff().up;
  27. // _isstemp.str(str);
  28. // for (double temp; _isstemp>>temp;)
  29. // x.push_back(temp);
  30. //чтение нач. коорд. преобразование для ginac согласно 5.15.2 Expression input
  31. for (unsigned i=0;i<x.size();i++)
  32. {
  33. str=str_sm + "_"+ std::to_string(i)+ "_";
  34. var_name.push_back(GiNaC::symbol(str));
  35. table[str]=var_name[i];
  36. }
  37. // var_name.push_back(GiNaC::symbol("a"));
  38. //table["a"]=var_name.back();//проверку на наличия сделать здесь
  39. //for (std::vector<std::string>::iterator it = uprav_par.begin(); it != uprav_par.end(); it++)
  40. // table[*it]=GiNaC::symbol(*it);
  41. table["xy"]=GiNaC::symbol("xy");//в функцую + проверять есть circle3
  42. table["xz"]=GiNaC::symbol("xz");
  43. table["yx"]=GiNaC::symbol("yx");
  44. table["yz"]=GiNaC::symbol("yz");
  45. table["zx"]=GiNaC::symbol("zx");
  46. table["zy"]=GiNaC::symbol("zy");
  47. for (unsigned i=0;i<uprav_par_naz.size();i++)
  48. {
  49. uprav_par_ex.push_back(GiNaC::symbol(uprav_par_naz[i]));
  50. table[uprav_par_naz[i]]=uprav_par_ex[i];
  51. }
  52. // GiNaC::ex cc=GiNaC::symbol("a");
  53. // table["a"]=cc;
  54. GiNaC::parser reader(table, true); //если есть различие между переменными в уравнениях и объявлении, то true выдает ошибку
  55. equations=sod->get_par_diff().eqn;
  56. equations_svob=sod->get_par_diff().eqn_svob;
  57. est_uprav=false; //по умолчанию считаем
  58. for (auto it = equations.begin(); it != equations.end(); it++)
  59. {
  60. expr_in.push_back(reader(*it));
  61. for (size_t j=0; j<uprav_par_ex.size(); j++)
  62. est_uprav=est_uprav || GiNaC::has(expr_in.back(),uprav_par_ex[j]);
  63. }
  64. equations.insert(equations.end(), equations_svob.begin(), equations_svob.end());
  65. for (auto it = equations_svob.begin(); it != equations_svob.end(); it++)
  66. {
  67. expr_svob.push_back(reader(*it));
  68. for (size_t j=0; j<uprav_par_ex.size(); j++)
  69. est_uprav=est_uprav || GiNaC::has(expr_svob.back(),uprav_par_ex[j]);
  70. }
  71. if (x.size()-expr_in.size() != 1 && (expr_svob.size() == 0 || expr_svob.size() % (x.size()-expr_in.size()-1) != 0))
  72. {
  73. _isstemp.str("");
  74. _isstemp << "Неправильное число данных. Число начальных координат " << x.size() << " число постоянных уравнений " << expr_in.size() << " число уравнений для степени свободы больше одной " << expr_svob.size();
  75. throw std::invalid_argument(_isstemp.str());
  76. }
  77. n_svob=sod->get_par_diff().nsv;
  78. N_var=sod->get_par_diff().nv;
  79. N_ex=expr_in.size(); //число уравнений
  80. if (!sod->get_par_diff().npp)
  81. det_comp();
  82. else
  83. {
  84. std::ifstream _cppsdu("test.cpp");
  85. if (_cppsdu.fail())
  86. throw std::invalid_argument("Файл сду cpp отсутствует");
  87. }
  88. }
  89. void sys_diff_equation::det_comp()
  90. {
  91. std::ofstream ofst{"test.cpp"}; // текстовая форма сду
  92. ofst << "#include <cmath> " << std::endl;
  93. ofst << "#include \"virteqn.h\"" << '\n' << std::endl;
  94. ofst << "class funcinst : public diffeqninterface {" << std::endl << "public:" <<std::endl;
  95. ofst << "virtual std::vector<syseqn> funcv();" << std::endl;
  96. ofst << "virtual state_type get_upr_par()";
  97. if (!est_uprav)
  98. ofst << " {return {};}" << std::endl;
  99. else
  100. ofst << ";" << std::endl;
  101. ofst << "private:" << std::endl;
  102. for (auto it = uprav_par_naz.begin(); it != uprav_par_naz.end(); it++)
  103. ofst << "static double "<< *it <<";" << std::endl;
  104. //ofst << "virtual void set_par_a(double &_a) {a=_a;}" <<std::endl;
  105. GiNaC::ex temp_1=expr_in.back();
  106. //bool check_t=false;
  107. //плохой код, доработать
  108. // for (int i=0; i<N_ex;i++)
  109. // if (expr_in[i].diff(var_name.back()) != 0)
  110. // check_t=true;
  111. // if (check_t) N_var++;
  112. for (int k=0; k < n_svob;k++)
  113. ofst << "static void func_" << k+1 << "(const state_type &, state_type &, const double);" << std::endl;
  114. ofst << "};" << std::endl << std::endl;
  115. ofst << "extern \"C\" diffeqninterface* create() {" << "return new funcinst;" << "}"<<std::endl;
  116. ofst << "extern \"C\" void destroy(diffeqninterface *p) {" << "delete p;" << "}"<<std::endl <<std::endl;
  117. for (auto it = uprav_par_naz.begin(); it != uprav_par_naz.end() && est_uprav; it++)
  118. ofst << "double funcinst::"<< *it << "=0;" << std::endl;
  119. if (est_uprav)
  120. {
  121. ofst << "state_type funcinst::get_upr_par() {" << std::endl << "state_type _upr_par;" << std::endl;
  122. for (auto it = uprav_par_naz.begin(); it != uprav_par_naz.end() && est_uprav; it++)
  123. ofst << "_upr_par.push_back("<< *it << ");" << std::endl;
  124. ofst << "return _upr_par;" << std::endl << "}" << std::endl << std::endl;
  125. }
  126. ofst << "std::vector<syseqn> funcinst::funcv() {" << std::endl;
  127. ofst << "std::vector<syseqn> temp_;" << std::endl;
  128. for (int k=0; k < n_svob;k++)
  129. ofst << "temp_.push_back(func_" << k+1 << ");" << std::endl;
  130. ofst << "return temp_;" << std::endl << "}" << std::endl << std::endl;
  131. for (int k=0; k < n_svob;k++)
  132. {
  133. if (n_svob != 1)
  134. expr_in.push_back(expr_svob[k]);
  135. GiNaC::lst biglst, lst_svb;
  136. for (int j=0; j<N_var-1;j++)
  137. for (int i=0; i<N_var-1;i++)
  138. biglst.append(expr_in[i].diff(var_name[j]));
  139. //столбец свободных членов со знаком "-"
  140. for (int i=0; i<N_var-1;i++)
  141. lst_svb.append(-expr_in[i].diff(var_name[N_var-1]));
  142. //матрица имеет N_var*N_var, но счет начинается с 0, поэтому N_var-1
  143. expr_in.pop_back();
  144. GiNaC::matrix A_mat(N_var-1, N_var-1, biglst);
  145. GiNaC::matrix b_mat(1, N_var-1, lst_svb), b_temp(1, N_var-1, lst_svb);
  146. ofst << "void funcinst::func_"<<k+1<<"(const state_type &x, state_type &f, const double t)" << std::endl;
  147. ofst << "{" << std::endl;
  148. for (int count=0; count<N_var-1; ++count)
  149. {
  150. ofst << "f[" << count << "] = (";
  151. //ginac нету функции менять столбец, удаление столбца для четных значений приводит к неправильному знаку
  152. for (int i=0; i<N_var-1;i++)
  153. {
  154. b_temp(0,i)=A_mat(count,i);
  155. A_mat(count,i)=b_mat(0,i);
  156. }
  157. std::stringstream prav_ch; //правая часть
  158. GiNaC::determinant(A_mat).print(GiNaC::print_csrc_double(prav_ch));
  159. ofst << vnut_pred(prav_ch).rdbuf() << ");" << std::endl;
  160. for (int i=0; i<N_var-1;i++)
  161. A_mat(count,i)=b_temp(0,i);
  162. }
  163. ofst << "f[" << N_var-1 << "] = (";
  164. std::stringstream prav_ch; //правая часть
  165. GiNaC::determinant(A_mat).print(GiNaC::print_csrc_double(prav_ch));
  166. ofst << vnut_pred(prav_ch).rdbuf() << ");" << std::endl << "}" << std::endl;
  167. }
  168. if (n_svob == 1)
  169. expr_in.push_back(temp_1);
  170. ofst.close();
  171. }
  172. ginacdata sys_diff_equation::get_ginacdata()
  173. {
  174. ginacdata A;
  175. A.vpnum=vp_num;
  176. A.cnum=circ_num;
  177. A.ccord=circ_coord;
  178. return A;
  179. }
  180. void sys_diff_equation::info_eqn_octave(std::stringstream &error_stream, std::stringstream &eqn_stream)
  181. {
  182. eqn_stream<<"#{\n";
  183. std::regex regSn(":"), regz("\\*|\\^|\\/");
  184. for (size_t i=0; i<expr_in.size(); i++)
  185. {
  186. std::stringstream strsm,strsm_;
  187. expr_in[i].print(GiNaC::print_dflt(strsm_));
  188. strsm << vnut_pred(strsm_).rdbuf();
  189. std::string str(strsm.str());
  190. const char coord_name[]="xyz";
  191. std::size_t pos_l = str.find('[', 0)+1;//следующий символ цифра
  192. std::size_t pos_r = str.find(']', 0)+1;
  193. while( pos_l !=std::string::npos+1)
  194. {
  195. unsigned _numb=std::stoi(str.substr(pos_l, pos_r-pos_l));
  196. std::string _str_r;
  197. _str_r=coord_name[_numb % dim_line];
  198. _str_r+="(:," + std::to_string(_numb/dim_line+1) + ")";
  199. str.replace(pos_l-2, pos_r-(pos_l-2), _str_r);//начиная с x[ т.е -2
  200. pos_l = str.find('[', pos_l+1)+1;
  201. pos_r = str.find(']', pos_r+1)+1;
  202. }
  203. str=regex_replace(str, regz, ".$&");
  204. eqn_stream << "eqn(:," << i+1 << ")=" << str << ";" << std::endl;
  205. error_stream << "error_eqn(:," << i+1 << ")=" << str << " - (" <<regex_replace(str, regSn, "1") << ");" <<std::endl;
  206. }
  207. if (uprav_par_naz.size() == 0)
  208. {
  209. eqn_stream << "#}\n";
  210. error_stream << "max_error=max(abs(error_eqn));\n\n";
  211. }
  212. else
  213. error_stream << "max_error=max(abs(error_eqn));\n#}\n\n";
  214. }
  215. std::stringstream &sys_diff_equation::vnut_pred(std::stringstream &prav_ch)
  216. {
  217. std::string input=prav_ch.str();
  218. std::regex left_s("_([[:d:]]+)"), right_s("([[:d:]]+)_");
  219. input=regex_replace(input,left_s,"[$1");
  220. input=regex_replace(input,right_s,"$1]");
  221. prav_ch.str("");
  222. prav_ch << input;
  223. return prav_ch;
  224. }