SHARE
TWEET

david

shadeyourself Jun 5th, 2020 911 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. double fST(double s, double t){
  5.     return 1/sqrt(pow(s, 4) + pow(t, 4));
  6. }
  7.  
  8. double integrate_dt(double s, double a, double b, double n){
  9.     double h = (b-a)/n;
  10.     double sum = 0;
  11.  
  12.     for (int i = 1; i <= n / 2; i++)
  13.         sum += 4 * fST(s, a + (2 * i - 1) * h);
  14.  
  15.     for (int i = 1; i < n / 2; i++)
  16.         sum += 2 * fST(s, a + 2 * i * h);
  17.  
  18.     return (fST(s, a) + sum + fST(s, b)) * h / 3;
  19. }
  20. double get_fS(double s, double a, double b, long n, double eps) {
  21.     double x = integrate_dt(s, a, b, n);
  22.     double y = integrate_dt(s, a, b, 2 * n);
  23.  
  24.     if (fabs(x - y) < eps)
  25.         return y;
  26.  
  27.     return get_fS(s, a, (a+b)/2, n, eps) + get_fS(s, (a + b) / 2, b, n, eps);
  28. }
  29. double fS(double s, double eps){
  30.     return get_fS(s, 0, 1, 2, eps);
  31. }
  32. //check fS value when s is some const here
  33. //сheking succseed
  34.  
  35. //now we have some f(s) func
  36. //the problem is now to find x for integral f(s) ds from 1 to x = alpha
  37.  
  38. double integrate(double a, double b, long n, double eps){
  39.     double h = (b-a)/n;
  40.     double sum = 0;
  41.  
  42.     for (int i = 1; i <= n / 2; i++)
  43.         sum += 4 * fS(a + (2 * i - 1) * h, eps);
  44.  
  45.     for (int i = 1; i < n / 2; i++)
  46.         sum += 2 * fS(a + 2 * i * h, eps);
  47.  
  48.     return (fS(a, eps) + sum + fS(b, eps)) * h / 3;
  49. }
  50. double getIntegral(double a, double b, long n, double eps) {
  51.     double x = integrate(a, b, n, eps);
  52.     double y = integrate(a, b, 2 * n, eps);
  53.  
  54.     if (fabs(x - y) < eps)
  55.         return y;
  56.     //printf("!\n");
  57.     return getIntegral(a, (a+b)/2, n, eps) + getIntegral((a + b) / 2, b, n, eps);
  58. }
  59. double f(double x, double alpha, double eps){
  60.     return (getIntegral(1, x, 2, eps) - alpha);
  61. }
  62. //check f + alpha value
  63. //checking succeed
  64.  
  65. //f is now a func we need an x for
  66. //so lets make a dev for it and then use the newton method
  67.  
  68. double df(double x, double alpha, double eps){
  69.     return (f(x+eps, alpha, eps) - f(x-eps, alpha, eps))/(2*eps);
  70. }
  71. //check df value
  72.  
  73. double newtone_method(double x0, double alpha, double eps) {
  74.     double x;
  75.     do {
  76.         x = x0;
  77.         x0 -= f(x0, alpha, eps) / df(x0, alpha, eps);
  78.     } while (fabs(x - x0) > eps);
  79.     return x0;
  80. }
  81.  
  82. int main(){
  83.     double alpha;
  84.     double eps;
  85.     printf("input alpha: ");
  86.     scanf("%lf", &alpha);
  87.     if(alpha >=1 || alpha <0){
  88.         printf("invalid alpha: integral cant be 1 or more or less than 0\n");
  89.         return -1;
  90.     }
  91.     printf("input eps: ");
  92.     scanf("%lf", &eps);
  93.     if(eps > 1 || eps < 1e-8){
  94.         printf("invalid eps\n");
  95.         return -1;
  96.     }
  97.     double x = newtone_method(1, alpha, eps);
  98.     printf("x = %lf\nf(x) = %lf", x, f(x+eps, alpha, eps));
  99.  
  100.     //checking fS
  101.     //printf("%lf", fS(1, eps));
  102.     //UPD succeed
  103.  
  104.     //checking f
  105.     //printf("%lf\n", f(1+eps, alpha, eps));
  106.     //UPD succeed
  107.     //cheking df
  108.     //printf("%lf\n", df(2, alpha, eps));
  109.     //printf("%lf\n", fS(2, eps));
  110.     //UPD succeed
  111.     return 0;
  112. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top