Dichotomy_LGR.cpp 1.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  1. // 输入: Pos:1 * 2,两端点
  2. // n:阶数
  3. // 输出: tn:精确的零点
  4. // 方法: 二分法(LGR)
  5. // 编写: 李兆亭,2019 / 11 / 04
  6. #include "GPM.h"
  7. int Dichotomy_LGR(int* n, vec* Pos, double* tn) {
  8. double tol = 1e-12;
  9. int Iscontinue = 1;
  10. double fx1, fx2, fx;
  11. double fy1, fy2, fy;
  12. double fav1, fav2, fav;
  13. double x = (*Pos)(0);
  14. double y = (*Pos)(1);
  15. double av = (x + y) / 2;
  16. int n1 = (*n) + 1;
  17. Legendre(&n1, &x, &fx1);
  18. Legendre(n, &x, &fx2);
  19. fx = fx1 + fx2;
  20. // cout << fx << endl;
  21. Legendre(&n1, &y, &fy1);
  22. Legendre(n, &y, &fy2);
  23. fy = fy1 + fy2;
  24. Legendre(&n1, &av, &fav1);
  25. Legendre(n, &av, &fav2);
  26. fav = fav1 + fav2;
  27. while (Iscontinue == 1) {
  28. if (abs(fx) <= tol) {
  29. *tn = x;
  30. Iscontinue = 0;
  31. }
  32. else if (abs(fy) <= tol) {
  33. *tn = y;
  34. Iscontinue = 0;
  35. }
  36. else if (abs(fav) <= tol) {
  37. *tn = av;
  38. Iscontinue = 0;
  39. }
  40. else {
  41. if (fx*fav < 0) {
  42. x = x; y = av; av = (x + y) / 2;
  43. Legendre(&n1, &y, &fy1);
  44. Legendre(n, &y, &fy2);
  45. fy = fy1 + fy2;
  46. Legendre(&n1, &av, &fav1);
  47. Legendre(n, &av, &fav2);
  48. fav = fav1 + fav2;
  49. }
  50. else if (fav*fy < 0) {
  51. x = av; y = y; av = (x + y) / 2;
  52. Legendre(&n1, &x, &fx1);
  53. Legendre(n, &x, &fx2);
  54. fx = fx1 + fx2;
  55. Legendre(&n1, &av, &fav1);
  56. Legendre(n, &av, &fav2);
  57. fav = fav1 + fav2;
  58. }
  59. }
  60. }
  61. return 1;
  62. }