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

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