Решение систем нелинейных уравнений https://www.mapleprimes.com/users/one%20man/posts?page=1
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

sysdiffeqn.cpp 9.9KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261
  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. sod->est_upr_par(est_uprav);
  72. if (x.size()-expr_in.size() != 1 && (expr_svob.size() == 0 || expr_svob.size() % (x.size()-expr_in.size()-1) != 0))
  73. {
  74. _isstemp.str("");
  75. _isstemp << "Неправильное число данных. Число начальных координат " << x.size() << " число постоянных уравнений " << expr_in.size() << " число уравнений для степени свободы больше одной " << expr_svob.size();
  76. throw std::invalid_argument(_isstemp.str());
  77. }
  78. n_svob=sod->get_par_diff().nsv;
  79. N_var=sod->get_par_diff().nv;
  80. N_ex=expr_in.size(); //число уравнений
  81. if (!sod->get_par_diff().npp)
  82. det_comp();
  83. else
  84. {
  85. std::ifstream _cppsdu("test.cpp");
  86. if (_cppsdu.fail())
  87. throw std::invalid_argument("Файл сду cpp отсутствует");
  88. }
  89. }
  90. void sys_diff_equation::det_comp()
  91. {
  92. std::ofstream ofst{"test.cpp"}; // текстовая форма сду
  93. if (uprav_par_naz.size()>0 && !est_uprav)
  94. {
  95. ofst <<"//Не используется управляемый параметр" << std::endl;
  96. // _isstemp.str("");
  97. // _isstemp << "Не используется управляемый параметр";
  98. // throw std::invalid_argument(_isstemp.str());
  99. }
  100. ofst << "#include <cmath> " << std::endl;
  101. ofst << "#include \"virteqn.h\"" << '\n' << std::endl;
  102. ofst << "class funcinst : public diffeqninterface {" << std::endl << "public:" <<std::endl;
  103. ofst << "virtual std::vector<syseqn> funcv();" << std::endl;
  104. ofst << "virtual state_type get_upr_par()";
  105. // if (!est_uprav)
  106. // ofst << " {return {};}" << std::endl;
  107. //else
  108. ofst << ";" << std::endl;
  109. ofst << "private:" << std::endl;
  110. for (auto it = uprav_par_naz.begin(); it != uprav_par_naz.end(); it++)
  111. ofst << "static double "<< *it <<";" << std::endl;
  112. //ofst << "virtual void set_par_a(double &_a) {a=_a;}" <<std::endl;
  113. GiNaC::ex temp_1=expr_in.back();
  114. //bool check_t=false;
  115. //плохой код, доработать
  116. // for (int i=0; i<N_ex;i++)
  117. // if (expr_in[i].diff(var_name.back()) != 0)
  118. // check_t=true;
  119. // if (check_t) N_var++;
  120. for (int k=0; k < n_svob;k++)
  121. ofst << "static void func_" << k+1 << "(const state_type &, state_type &, const double);" << std::endl;
  122. ofst << "};" << std::endl << std::endl;
  123. ofst << "extern \"C\" diffeqninterface* create() {" << "return new funcinst;" << "}"<<std::endl;
  124. ofst << "extern \"C\" void destroy(diffeqninterface *p) {" << "delete p;" << "}"<<std::endl <<std::endl;
  125. for (auto it = uprav_par_naz.begin(); it != uprav_par_naz.end(); it++)
  126. ofst << "double funcinst::"<< *it << "=0;" << std::endl;
  127. // if (est_uprav)
  128. // {
  129. ofst << "state_type funcinst::get_upr_par() {" << std::endl << "state_type _upr_par;" << std::endl;
  130. for (auto it = uprav_par_naz.begin(); it != uprav_par_naz.end(); it++)
  131. ofst << "_upr_par.push_back("<< *it << ");" << std::endl;
  132. ofst << "return _upr_par;" << std::endl << "}" << std::endl << std::endl;
  133. // }
  134. ofst << "std::vector<syseqn> funcinst::funcv() {" << std::endl;
  135. ofst << "std::vector<syseqn> temp_;" << std::endl;
  136. for (int k=0; k < n_svob;k++)
  137. ofst << "temp_.push_back(func_" << k+1 << ");" << std::endl;
  138. ofst << "return temp_;" << std::endl << "}" << std::endl << std::endl;
  139. for (int k=0; k < n_svob;k++)
  140. {
  141. if (n_svob != 1)
  142. expr_in.push_back(expr_svob[k]);
  143. GiNaC::lst biglst, lst_svb;
  144. for (int j=0; j<N_var-1;j++)
  145. for (int i=0; i<N_var-1;i++)
  146. biglst.append(expr_in[i].diff(var_name[j]));
  147. //столбец свободных членов со знаком "-"
  148. for (int i=0; i<N_var-1;i++)
  149. lst_svb.append(-expr_in[i].diff(var_name[N_var-1]));
  150. //матрица имеет N_var*N_var, но счет начинается с 0, поэтому N_var-1
  151. expr_in.pop_back();
  152. GiNaC::matrix A_mat(N_var-1, N_var-1, biglst);
  153. GiNaC::matrix b_mat(1, N_var-1, lst_svb), b_temp(1, N_var-1, lst_svb);
  154. ofst << "void funcinst::func_"<<k+1<<"(const state_type &x, state_type &f, const double t)" << std::endl;
  155. ofst << "{" << std::endl;
  156. for (int count=0; count<N_var-1; ++count)
  157. {
  158. ofst << "f[" << count << "] = (";
  159. //ginac нету функции менять столбец, удаление столбца для четных значений приводит к неправильному знаку
  160. for (int i=0; i<N_var-1;i++)
  161. {
  162. b_temp(0,i)=A_mat(count,i);
  163. A_mat(count,i)=b_mat(0,i);
  164. }
  165. std::stringstream prav_ch; //правая часть
  166. GiNaC::determinant(A_mat).print(GiNaC::print_csrc_double(prav_ch));
  167. ofst << vnut_pred(prav_ch).rdbuf() << ");" << std::endl;
  168. for (int i=0; i<N_var-1;i++)
  169. A_mat(count,i)=b_temp(0,i);
  170. }
  171. ofst << "f[" << N_var-1 << "] = (";
  172. std::stringstream prav_ch; //правая часть
  173. GiNaC::determinant(A_mat).print(GiNaC::print_csrc_double(prav_ch));
  174. ofst << vnut_pred(prav_ch).rdbuf() << ");" << std::endl << "}" << std::endl;
  175. }
  176. if (n_svob == 1)
  177. expr_in.push_back(temp_1);
  178. ofst.close();
  179. }
  180. ginacdata sys_diff_equation::get_ginacdata()
  181. {
  182. ginacdata A;
  183. A.vpnum=vp_num;
  184. A.cnum=circ_num;
  185. A.ccord=circ_coord;
  186. return A;
  187. }
  188. void sys_diff_equation::info_eqn_octave(std::stringstream &error_stream, std::stringstream &eqn_stream)
  189. {
  190. eqn_stream<<"#{\n";
  191. std::regex regSn(":"), regz("\\*|\\^|\\/");
  192. for (size_t i=0; i<expr_in.size(); i++)
  193. {
  194. std::stringstream strsm,strsm_;
  195. expr_in[i].print(GiNaC::print_dflt(strsm_));
  196. strsm << vnut_pred(strsm_).rdbuf();
  197. std::string str(strsm.str());
  198. const char coord_name[]="xyz";
  199. std::size_t pos_l = str.find('[', 0)+1;//следующий символ цифра
  200. std::size_t pos_r = str.find(']', 0)+1;
  201. while( pos_l !=std::string::npos+1)
  202. {
  203. unsigned _numb=std::stoi(str.substr(pos_l, pos_r-pos_l));
  204. std::string _str_r;
  205. _str_r=coord_name[_numb % dim_line];
  206. _str_r+="(:," + std::to_string(_numb/dim_line+1) + ")";
  207. str.replace(pos_l-2, pos_r-(pos_l-2), _str_r);//начиная с x[ т.е -2
  208. pos_l = str.find('[', pos_l+1)+1;
  209. pos_r = str.find(']', pos_r+1)+1;
  210. }
  211. str=regex_replace(str, regz, ".$&");
  212. eqn_stream << "eqn(:," << i+1 << ")=" << str << ";" << std::endl;
  213. error_stream << "error_eqn(:," << i+1 << ")=" << str << " - (" <<regex_replace(str, regSn, "1") << ");" <<std::endl;
  214. }
  215. if (uprav_par_naz.size() == 0)
  216. {
  217. eqn_stream << "#}\n";
  218. error_stream << "max_error=max(abs(error_eqn));\n\n";
  219. }
  220. else
  221. error_stream << "max_error=max(abs(error_eqn));\n#}\n\n";
  222. }
  223. std::stringstream &sys_diff_equation::vnut_pred(std::stringstream &prav_ch)
  224. {
  225. std::string input=prav_ch.str();
  226. std::regex left_s("_([[:d:]]+)"), right_s("([[:d:]]+)_");
  227. input=regex_replace(input,left_s,"[$1");
  228. input=regex_replace(input,right_s,"$1]");
  229. prav_ch.str("");
  230. prav_ch << input;
  231. return prav_ch;
  232. }