Решение систем нелинейных уравнений 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.

octaveplg.cpp 6.4KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. #include "octaveplg.h"
  2. octaveplg::octaveplg(sys_osn_dan &sod_, sys_diff_equation &sde_, std::string &name_, bool(_bench_m)): sod(&sod_), sde(&sde_), octout(name_), bench_m(_bench_m)
  3. {
  4. x=sod->get_par_diff().x0;
  5. n_svob=sod->get_par_diff().nsv;
  6. equations=sod->get_par_diff().eqn;
  7. dim_line=sod->get_par_diff().dln;
  8. est_upr=sod->get_par_diff().up.size();
  9. vp_num=sde->get_ginacdata().vpnum;
  10. circ_num=sde->get_ginacdata().cnum;
  11. circ_coord=sde->get_ginacdata().ccord;
  12. sde->info_eqn_octave(error_stream, eqn_stream);
  13. }
  14. void octaveplg::octave_m()
  15. {
  16. const char *coord_nam="xyz";
  17. out << "#{\n";
  18. for (auto i: equations)
  19. out << i << "\n";
  20. for (auto i: x)
  21. out << (i) << " ";
  22. out << "\n#}\n\n" << "clear all;\n\n";
  23. if (bench_m)
  24. out << "#b_var=logspace(,);\n" << "for iter=1:columns(b_var)\n" <<"clear sol error_eqn;\n";
  25. out << "tic();\n";
  26. if (est_upr != 0)
  27. out << "[sol dt upr]=";
  28. else
  29. out << "[sol dt]=";
  30. if (bench_m)
  31. out << "draghilev_oct(-1,b_var(iter));\n";
  32. else
  33. out << "draghilev_oct;\n";
  34. out<< "draghilev_time=toc();\n" << "num=rows(sol);\n" << "N_col=columns(sol);\n";
  35. if (n_svob != 1)
  36. {
  37. out << "if (mod(num, " << n_svob <<") != 0)\n num=num-mod(num," << n_svob << ");\n" << "endif\n\n";
  38. out << "sol=sol(1:num,:);\n";
  39. for (int i=0;i<dim_line;i++)
  40. out << coord_nam[i] << " = (sol(1:" << n_svob << ":end, " << i+1 << ":" << dim_line << ":N_col) + sol(2:" << n_svob << ":end, " << i+1 << ":" << dim_line << ":N_col))/" << n_svob << ";\n";
  41. }
  42. else
  43. for (int i=0;i<dim_line;i++)
  44. out << coord_nam[i] << " = sol(:, " << i+1 << ":" << dim_line << ":N_col);\n";
  45. out << "\n"<< eqn_stream.rdbuf();
  46. out << "\n"<< error_stream.rdbuf();
  47. if (bench_m)
  48. out << "err_bench(iter,:)=max_error;\n" << "time_bench(iter,:)=draghilev_time;\n" <<"endfor\n\n" << "loglog(time_bench, err_bench, \"marker\", \"o\", \"linestyle\", \"none\");\n"<<"xlabel(\"Время, с\");\n"<<"ylabel(\"Максимальная абсолютная погрешность по модулю\");\n"<<"set (gca, \"ygrid\", \"on\");";
  49. else
  50. write_graphic_tom();
  51. octout<<out.str();
  52. octout.close();
  53. }
  54. void octaveplg::write_graphic_tom()//to m-script
  55. {
  56. std::vector<const char*> colour_prt ={ " \"color\", \"r\", \"marker\", \"none\",\n",
  57. " \"color\", \"g\", \"marker\", \"none\",\n",
  58. " \"b\", \"marker\", \"o\",\n",
  59. " \"b\", \"marker\", \"+\",\n",
  60. "\"color\", \"blue\",",
  61. "\"color\", \"cyan\"," };
  62. std::vector<const char*> prt_point ={ "1:m,", "m:k,", "k," }; //параметры для точки
  63. int N_var=x.size();
  64. int n_point=N_var/dim_line;
  65. int nln=n_point*prt_point.size()+circ_num.size()+1; //номер 1-ой линии
  66. const char *coord_nam="xyz";//повтор, нужно убрать
  67. for (unsigned i=0;i<circ_num.size();i++) //1 точка 1 параметр для окружности
  68. {
  69. for (int k=0;k<dim_line;k++)
  70. out << coord_nam[k] << i << "=" <<circ_coord[dim_line*i+k]<< ";\t";
  71. out << "\n";
  72. }
  73. out << "h = max(mean(abs(diff(sol))))/10;\n";
  74. out << "limits = [";
  75. for (int i=0;i<dim_line;i++)
  76. out << "min(" << coord_nam[i] << "(:)), max(" << coord_nam[i] << "(:))+h, ";
  77. out.seekp(-2, std::ios::cur);
  78. out << "];\n" << "limits = [";
  79. for (int i=0;i<dim_line;i++)
  80. {
  81. out << "min([limits(" << 2*i+1 << ") ";
  82. for (unsigned j=0;j<circ_num.size();j++)
  83. out << coord_nam[i] << j << " ";
  84. out.seekp(-1, std::ios::cur);
  85. out <<"]), max([limits(" << 2*i+2 << ") ";
  86. for (unsigned j=0;j<circ_num.size();j++)
  87. out << coord_nam[i] << j << " ";
  88. out.seekp(-1, std::ios::cur);
  89. out << "]) ,";
  90. }
  91. out.seekp(-2, std::ios::cur);
  92. out << "];\n\n";
  93. out << "hl = plot";
  94. if (dim_line ==3 )
  95. out << 3; //plot3 для 3d
  96. out <<"(\n";
  97. for (int i=0;i<n_point;i++)
  98. {
  99. for (unsigned j=0;j<prt_point.size();j++) //1 точка 3 параметра для линии
  100. {
  101. for (int k=0;k<dim_line;k++)
  102. out << "\t " << coord_nam[k] << "(1," << i+1 << "),";
  103. out << colour_prt[j];
  104. }
  105. }
  106. for (unsigned i=0;i<circ_num.size();i++) //1 точка 1 параметр для окружности
  107. {
  108. for (int k=0;k<dim_line;k++)
  109. out << "\t " << coord_nam[k] << i << ",";
  110. out << colour_prt[3];
  111. }
  112. for (auto i = vp_num.begin(); i != vp_num.end(); i+=2)//линия
  113. {
  114. for (int j=0;j<dim_line;j++)
  115. out << "\t ["<< coord_nam[j] << "(1," << *i << ") " << coord_nam[j] << "(1," << *(i+1) << ")], ";
  116. out << colour_prt[4] << "\n";
  117. }
  118. for (unsigned i=0;i<circ_num.size();i++)
  119. {
  120. for (int j=0;j<dim_line;j++)
  121. out << "\t ["<< coord_nam[j] << "(1," << circ_num[i] << ") " << coord_nam[j] << i << "], ";
  122. out << colour_prt[5] << "\n";
  123. }
  124. out.seekp(-2, std::ios::cur); //интересно, -2 отнимает два последних элемента c конца
  125. out << ");\n\n";
  126. if (n_svob != 1)
  127. out << "num=num/" << n_svob << ";\n";
  128. out << "dn = round (num/10);\n" << "#g_step=1;\n" << "g_step=round (num/120);\n" << "axis (\"equal\",1.001*limits);\n" << "axis (limits);\n" << "#legend({[";
  129. for (auto i: equations)
  130. out << "'" << i << "' char(10) ";
  131. out << "]}, 'location', 'northeast');\n" << "#legend boxoff\n\n";
  132. out << "for n = 1:g_step:(num+dn);\n" << "\t" << "m = n - dn;\n" << "\t" << "m = max ([m, 1]);\n" << "\t" << "k = min ([n, num]);\n";
  133. for (int i=0;i<n_point;i++)
  134. for (unsigned j=0;j<prt_point.size();j++)
  135. {
  136. out << "\t" << "set (hl(" << prt_point.size()*i+j+1 << "), ";
  137. for (int k=0;k<dim_line;k++)
  138. out << "\""<< coord_nam[k] <<"data\", " << coord_nam[k] << "(" << prt_point[j] << i+1 << "), ";
  139. out.seekp(-2, std::ios::cur);
  140. out << ");\n";
  141. }
  142. for (auto i = vp_num.begin(); i != vp_num.end(); i+=2)
  143. {
  144. out << "\t" << "set (hl(" << nln << "), ";
  145. for (int j=0;j<dim_line;j++)
  146. out << "\"" << coord_nam[j] << "data\", [" << coord_nam[j] << "(k," << *i << ") " << coord_nam[j] << "(k," << *(i+1) << ")], ";
  147. out.seekp(-2, std::ios::cur);
  148. out << ");\n";
  149. nln++;
  150. }
  151. for (unsigned i=0;i<circ_num.size();i++)
  152. {
  153. out << "\t" << "set (hl(" << nln << "), ";
  154. for (int j=0;j<dim_line;j++)
  155. out << "\"" << coord_nam[j] << "data\", [" << coord_nam[j] << "(k," << circ_num[i] << ") " << coord_nam[j] << i << "], ";
  156. out.seekp(-2, std::ios::cur);
  157. out << ");\n";
  158. nln++;
  159. }
  160. out << "\t" << "drawnow();\n" << "endfor";
  161. }