自适应辛普森(算法简要 + 模板)
簡述:
對于普通的二次函數有f(x)=ax2+bx+cf(x) = ax ^ 2 + bx + cf(x)=ax2+bx+c,原函數F(x)=ax33+bx22+cxF(x) = \frac{a x^3}{3} + \frac{bx ^2}{2} + cxF(x)=3ax3?+2bx2?+cx。
有
?lrf(x)=F(r)?F(l)=a3(r?l)3+b2(r?l)2+c(r?l)=(r?l)6(2a(l2+r2+lr)+3b(l+r)+c)=(r?l)6((al2+bl+c)+(ar2+br+c)+4(a((l+r)2)2+b(l+r)2+c))=(r?l)6(f(l)+f(r)+4f((l+r)2))\lmoustache_{l} ^{r} f(x) = F(r) - F(l) \\ = \frac{a}{3}(r - l) ^3 + \frac{b}{2} (r - l) ^2 + c(r - l)\\ = \frac{(r - l)}{6} (2a(l ^ 2 + r ^2 + lr) + 3b(l + r) + c)\\ = \frac{(r - l)}{6} ((al ^2 + bl + c) +( a r ^ 2 + br + c) + 4(a(\frac{(l + r)}{2}) ^ 2 + b\frac{(l + r)}{2} + c))\\ = \frac{(r - l)}{6} (f(l) + f(r) + 4f(\frac{(l + r)}{2}))\\ ?lr?f(x)=F(r)?F(l)=3a?(r?l)3+2b?(r?l)2+c(r?l)=6(r?l)?(2a(l2+r2+lr)+3b(l+r)+c)=6(r?l)?((al2+bl+c)+(ar2+br+c)+4(a(2(l+r)?)2+b2(l+r)?+c))=6(r?l)?(f(l)+f(r)+4f(2(l+r)?))
這是普通的辛普森推導,自適應辛普森就是在普通的辛普森上加了一個評估步驟,我們通過把一個區間二分,如果這個區間整體求得的值與一半加一半單獨求得的值是較為接近的這個時候我們就可以認為我們已經求得正確答案。
P4525 【模板】自適應辛普森法1
/*Author : lifehappy */ #pragma GCC optimize(2) #pragma GCC optimize(3) #include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f; const double eps = 1e-6;double a, b, c, d, l, r;double f(double x) {return (c * x + d) / (a * x + b); }double sim(double l, double r) {return (r - l) * (f(l) + f(r) + 4.0 * f((l + r) / 2.0)) / 6.0; }double asr(double l, double r, double eps, double ans) {double mid = (l + r) / 2;double ansl = sim(l, mid), ansr = sim(mid, r);if(fabs(ansl + ansr - ans) < 15 * eps) return ansl + ansr + (ansl + ansr - ans) / 15.0;return asr(l, mid, eps / 2.0, ansl) + asr(mid, r, eps / 2.0, ansr); }int main() {// freopen("in.txt", "r", stdin); // freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);scanf("%lf %lf %lf %lf %lf %lf", &a, &b, &c, &d, &l, &r);printf("%.6f\n", asr(l, r, eps, sim(l, r)));return 0; }Ellipse
顯然有y=b×1?(xa)2y = b \times \sqrt{1 - (\frac{x}{a}) ^2}y=b×1?(ax?)2?,所以只要簡單的帶入求解即可。
/*Author : lifehappy */ #pragma GCC optimize(2) #pragma GCC optimize(3) #include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f; const double eps = 1e-7;const int N = 1e3 + 10, mod = 6;double a, b, l, r;double f(double x) {return 2 * b * (sqrt(1 - (x * x) / (a * a))); }double sim(double l, double r) {return (r - l) * (f(l) + f(r) + 4.0 * f((l + r) / 2.0)) / 6.0; }double asr(double l, double r, double eps, double ans) {double mid = (l + r) / 2;double ansl = sim(l, mid), ansr = sim(mid, r);if(fabs(ansl + ansr - ans) < eps) return ansl + ansr;return asr(l, mid, eps / 2.0, ansl) + asr(mid, r, eps / 2.0, ansr); }int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int T; scanf("%d", &T);while(T--) {scanf("%lf %lf %lf %lf", &a, &b, &l, &r);printf("%.3f\n", asr(l, r, 1e-5, sim(l, r)));}return 0; }The area
兩個函數積分求差,也就是分別算出兩個積分,然后相減即可。
對于二次函數的式子,這里直接用了拉格朗日插值來求了,沒去推導公式。
對于一次函數的式子,只要簡單求解一下即可。
/*Author : lifehappy */ #pragma GCC optimize(2) #pragma GCC optimize(3) #include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f; const double eps = 1e-6;double x[5], y[5];double f1(double X) {double ans = 0;for(int i = 1; i <= 3; i++) {double a = 1, b = 1;for(int j = 1; j <= 3; j++) {if(i == j) continue;a *= X - x[j];b *= x[i] - x[j];}ans += y[i] * a / b;}return ans; }double sim1(double l, double r) {return (r - l) * (f1(l) + f1(r) + f1((l + r) / 2.0) * 4.0) / 6.0; }double asr1(double l, double r, double eps, double ans) {double mid = (l + r) / 2.0;double ansl = sim1(l, mid), ansr = sim1(mid, r);if(fabs(ansl + ansr - ans) < 15.0 * eps) return ansl + ansr + (ansl + ansr - ans) / 25.0;return asr1(l, mid, eps / 2.0, ansl) + asr1(mid, r, eps / 2.0, ansr); }double f2(double X) {return (y[3] - y[2]) / (x[3] - x[2]) * (X - x[3]) + y[3]; }double sim2(double l, double r) {return (r - l) * (f2(l) + f2(r) + f2((l + r) / 2.0) * 4.0) / 6.0; }double asr2(double l, double r, double eps, double ans) {double mid = (l + r) / 2;double ansl = sim2(l, mid), ansr = sim2(mid, r);if(fabs(ansl + ansr - ans) < 15.0 * eps) return ansl + ansr + (ansl + ansr - ans) / 15.0;return asr2(l, mid, eps / 2.0, ansl) + asr2(mid, r, eps / 2.0, ansr); }int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int T; scanf("%d", &T);while(T--) {for(int i = 1; i <= 3; i++) {scanf("%lf %lf", &x[i], &y[i]);}double ans = asr1(x[2], x[3], eps, sim1(x[2], x[3])) - asr2(x[2], x[3], eps, sim2(x[2], x[3]));printf("%.2f\n", ans);}return 0; }總結
以上是生活随笔為你收集整理的自适应辛普森(算法简要 + 模板)的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 怎样瘦胸最快最有效
- 下一篇: 2019年ICPC银川区域赛 Easy