• API
• FAQ
• Tools
• Archive
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.
Top